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Abstract. The equations for laminar boundary layer flow over a general smooth 
surface in three dimensions are analyzed in a normal coordinate system. The invariance 
properties of these equations are found using the concept of subtensors. The boundary 
layer equations are not tensor equations but subtensor equations. Conditions for the 
Cartesian form of the equations are given and a criterion for no secondary flow is found 
in terms of the geodesics of the body surface. The displacement effect of the boundary 
layer is also discussed. 

Introduction. In recent times interest in three-dimensional boundary layer flows 
has grown considerably. A review of the subject has been given by Sears [1]. Most of 
the work cited by Sears is concerned with the solutions to particular problems. To the 
author’s knowledge the only general discussions of boundary layers in three-dimensional 
flows are contained in the work of Howarth [2], Moore [3], and Hayes [4]. However, 
there is also some work of C. C. Lin contained in Chap. 18 of [5], in which the appropriate 
equations are derived without much discussion. This latter work had been overlooked 
by the previously mentioned authors. The method presented by Lin proves to be useful 
in discussing the invariance properties of the boundary layer equations. (C. C. Lin has 
informed me that the material on general boundary layer equations [5] was supplied 
to the late Professor A. D. Michal not as new results, but as a presentation of some 
previous work by Levi-Civita.) 

In this paper it is first shown that Lin’s approach carries over to a more general class 
of coordinate systems. Certain invariance properties follow from examining the resulting 
equations. Now it is known that the two-dimensional boundary layer equations are not 
tensor equations and one would expect the same result in three-dimensions (In [4], 
Hayes seems to imply the contrary). However, in the language of [6], the boundary 
layer equations are subtensor equations which means that they are invariant with 
respect to certain types of coordinate transformations. This does not contradict the work 
of Lagerstrom and Kaplun [7] and [8], where the non-tensor character of the boundary 
layer equations is used to define an “optimum coordinate system”’’. A similar procedure 
could be developed for three-dimensional flows, however a more general coordinate 
system than that used here would have to be considered. 

From the subtensor form of the equations, it is easily shown that they reduce to 
the “Cartesian form” (as is true for any two-dimensional flow) for any surface whose 
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Gaussian curvature is zero if appropriate coordinates are used. Howarth [2], concludes 
that this happens only for planes and cylinders. The curvature effects that Howarth 
describes include the effects of a curvilinear coordinate system. (Moore [3] conjectured 
that the Cartesian form applies to any curved surface.) 

Using the streamline coordinates as in [4], it follows that a boundary layer flow is 
essentially two-dimensional if there is no secondary flow. A simple criterion for no 
secondary flow is given in terms of the geodesics of the surface. Finally the displacement 
effect of the boundary layer is discussed. This is done in a different way than by Moore 
[9], and some differences between the two- and three-dimensional cases are discussed. 

Only laminar flows of an incompressible, non-conducting fluid are discussed. Some 
of the conclusions should hold in more general cases, however. Also body forces have 
been neglected and difficulties due to separation effects are not discussed. 

The boundary layer equations. In applying the boundary layer concept to the 
flow of a viscous fluid we have some given surface to consider, for example a solid body 
or an interface between two fluids. This basic surface is used to define a coordinate 
system: it is one of the coordinate surfaces. A convenient class of coordinate systems 
is the normal coordinate systems, [6], which is defined as follows. A given one-parameter 
family of surfaces in a Riemannian N-space has a family of orthogonal trajectories 
such that, under very general conditions, only one trajectory passes through each point 
of space. The parameter of the family of surfaces is denoted by x’ and on one of the 
surfaces’ a coordinate system x’, --- , x’ is set up. We shall be concerned with 
a Euclidean 3-space although the derivation of the boundary layer equations would 
proceed in the same way for more general spaces. 

The convention is adopted that Greek suffixes have the range 1, 2 and small Latin 
suffixes the range 1, 2, 3. The given surface over which we consider the boundary layer 
flow is the one on which the coordinates x“ are defined and for this surface x° = 0. 
Normal coordinate systems have the property that the distance expression is given by 


ds’ = Jap ax" dx® +. Ja3(dx°)” 


that is, 


There is no unique way of choosing our normal coordinate system since all that is re- 
quired is that the given surface be one of the family. A simple special case is the geodesic 
normal coordinate system in which the family of surfaces is obtained by measuring off 
constant distances from a given surface along the geodesics which cut the surface ortho- 
gonally. In this special case 


733 _ i. 


This is the coordinate system used by Lin [5]. 

The surface x* = 0 is, in general, a Riemannian 2-space. In this subspace we can 
carry out tensor operations and these will be intimately related to the tensor operations 
of the parent 3-space. We use a comma to denote covariant differentiation in the 3-space 
and a semi-colon for the same operation in the subspace. A set of quantities 7’, which 
transform according to the law 


T! = T,(d2x°/dx'*) 
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under a transformation of the form 
x" = f “(°) 


13 3 
x = 2< 


(1) 


is called a subtensor. If a tensor 7’, is split up into two groups 7, , 7’; then for (1) T. 
is a subtensor and 7, is a subinvariant, and similarly for higher order tensors. Also 
the Christoffel symbols, in which one or more of the indices have the value 3, are sub- 
tensors. A more detailed discussion of substensors is given in [6]. 

The derivation of the boundary layer equations is straightforward but tedious. The 
procedure here is essentially the same as Lin’s for the special case of a geodesic normal 
coordinate system and therefore will not be discussed in detail. (Some errors were found 
in the details of [5] but these only affected terms which drop out in the boundary layer 
approximation.) Starting from the Navier-Stokes equations in tensor form 


(du;/dt) + wu,;.; = vg’ u; ja ~~ Bs (2) 
u’', = 0, 


where z is the ratio of pressure to the constant density, the following steps are necessary. 


The momentum equations are split up into two groups (as illustrated above with 7;); 


the covariant derivatives are expressed in terms of the “sub-covariant derivatives’’; a 


transformation 

a 3 1/2 4 f.1/2 

¢=a2/v U,=4u,/*"’; (3) 
is applied; all quantities, including the metric tensor and Christoffel symbols, are ex- 
panded in a power series in +'””; the equations for the lowest order terms then yield the 


boundary layer equations 


(du, At) + wuss + (U3/g$?)(du,/df) = —m,4 + (Ua ae’) / 933 
dxr/at = 0 (4) 
ur’, + Co? + (0U;/0¢)/gs = 0, 
where 
Ca = (0933 /dx*)/2g33° 
93: = Jas(z , xz’, 0) 
and the metric tensor for (4) is a,,3 Where 
Gop * Gaz, 2,0. 
Phe boundary layer equations (4) reduce to those given by Lin for g;; = 1. (Note that, 


in the procedure outlined above, nothing is implied about the higher approximations 
obtained from the series expansion. Consideration of these involves significant difficulties, 
[13] 

Under transformations of the form (1) it is easily seen that equations (4) are invariant 
since U, and g,, are invariant. Thus the boundary layer equagions have subtensor 
character. This is not too surprising since what destroys the tensor character of this 


approximation is the transformation (3) and this is not affected by (1). Also in the sub- 
tensor form it becomes obvious that the boundary layer equations reduce to “Cartesian 
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form” for any surface of zero curvature (flat 2-space) in which the coordinates x“ are 
Cartesian coordinates provided that a geodesic normal coordinate system is used, i.e. 
933 = 1. This conclusion differs from Howarth’s [2], because he restricts his analysis 
(which uses a geodesic normal system) to coordinate systems x“ which can be Cartesian 
only for planes and cylinders. Because of the subtensor character of (4) there is complete 
freedom in the choice of «“. However, the choice of Cartesian coordinates x“ on a de- 
velopable surface may not always be the most advantageous. For example, the form 
of the terms z., which are determined by the external flow, may be complicated by 
such a choice. It must be remembered that the conclusion concerning when the boundary 
layer equations reduce to Cartesian form applies only to within the approximations 
of the standard theory. For example, for the flow over on open-ended cylinder the 
equations are in Cartesian form (for appropriate x“) but the usual boundary layer 
approximations become invalid as the distance from the leading edge increases. Finally, 
it may be noted, even for flow over a plane the equations will have curvature terms 
appearing if x“ are not Cartesian coordinates. 

The boundary layer equations (4) can be written in the more conventional form in 
terms of the physical components of the velocity. Before doing this let us specialize to 
orthogonal coordinates x" and change the notation for the metric tensor 


g: = H;, 


ae he. (no summation) (5) 


9ss(x', x , 0) = hs = gaz , 
where h, is not a function of x*. Also denote the physical components of the velocity 
by u, v, and w. Thus 


’ 


u u,/h, 
9 = @e/Ry 
w= Us/Nsz., 


since for (4) the metric tensor components are the lower case h’s. In terms of the physical 


components, using x, y, and z as coordinates, Eqs. (4) become 

u, + uu,/h, + vu,/he + uvhy,/hyhe — vhee/hihe + wu./hy = —2,/hy + vu,,/h3 

v, + wv,/h, + w,/he + who,/hihe — why,/hihe + we./hg = —2,/he + w.,/hy , (6) 
y 


gr, =0 
[((hohgu), + (hihgv),|/Aihe + w./hz = 0, 


where the subscripts /, x, y, and z denote partial differentiation. These reduce to Howarth’s 
equations [2] for h; = 1, i.e. a geodesic normal coordinate system. Hayes [4], gives these 
equations (if compressibility is neglected in his equations) except that he, in some way, 
allows A, to depend on z. 

Finally, it can be remarked that, to examine the flow in the neighborhood of a 
stagnation point, as Howarth has done [10], it is only necessary to introduce a Riemannian 
coordinate system for x“ with origin at the stagnation point. 

Secondary flow. Since the boundary layer equations are invariant under trans- 
formations of the surface coordinates x“, it is natural to look for coordinates which, 
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under certain conditions, simplify the general equations. An ingenious choice of co- 
ordinates was made by Hayes [4]. For a steady external flow the coordinate z’ is chosen 
along the streamlines of the external flow evaluated on the surface and x’ is along the 
orthogonal trajectories of these streamlines. If U and V are the physical components 
of the external flow, evaluated on the surface, in the x and y (or x’ and 2°) directions, 


then for this choice of coordinates V = 0 and 
—a,/h, = UU,/h, , ” 
(7) 
—r,/h, = —U’°k, 
where 
k = h,,/h,hs 
i.e. k is the geodesic curvature of the streamlines. Consider steady flow in the boundary 
layer. The boundary conditions for v are 
v= 0, r—7 @ 
and, if the surface over which the flow is considered is a solid, non-spinning body, 
v = 0, z= 0. 


Now if & = 0, it is seen from (7) and the second equation of (6) that v = 0 is a solution. 
The term “secondary flow’ (also cross flow) is used to indicate that the streamlines in 
the boundary layer flow do not coincide with the external flow streamlines evaluated 
on the body. Also k = O has a simple geometric interpretation: the streamlines are 
geodesics of the surface. Thus, for steady flow over a non-spinning body, there is no 
secondary flow if the external flow streamlines are geodesics of the body surface. 

A simple example of the above conclusion is provided by the flow over a flat plate 
with an arbitrary leading edge placed in a uniform stream at no angle of attack; see 
Fig. 1. The external flow streamlines are just straight lines on the plate so there is no 


LEADING 
EDGE 


STREAMLINES 


Fia. 1. 
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secondary flow. (Moore [3], reaches this conclusion by finding a solution of the flow 
equations. He also indicates how difficulties arise if the leading edge does not have a 
continuously turning tangent.) A special case of this is the yawed flat plate. However 
consider yawed infinite cylinders, which many authors have done, [1]. It is easy to see, 
using the above criterion, that the yawed flat plate is the only case with no secondary 
flow since for any other cylinder the external flow streamlines cannot be geodesics. 
Since v = O in the streamline coordinates for k = 0, the boundary layer flow is 


essentially two-dimensional. Equations (6), written for a geodesic normal system, become 
uu,/h, + wu, = UU,/h, + m0.. , 
Tr, = Q, (8) 
(hou),/Ayhe + w 0. 


A new coordinate, NX, can be introduced, 


X=] hdx 
and then Eqs. (8) are seen to be in exactly the form of the boundary layer equations 
over a body of revolution. Therefore the same transformation that Mangler [11] has 
introduced will reduce equations (8) to the standard two-dimensional equations. 


Mangler’s transformation is 


é == | ea, #. @ 
n = hz, (9) 
n= 7 


w = how’ — hon’. 


Note that the coordinate transformation here is not of the type (1) and that the trans- 
formation to new velocity components (u’, w’) does not follow a tensor law. Hayes has 
indicated the possibility of transforming the no-secondary flow equations to two- 
dimensional form. However his transformation is made by choosing the form of the metric 
tensor in a suitable manner and it would be very difficult to find out exactly what the 
transformed coordinates are. From the discussion above it is seen that the transfor- 
mation is just that of Mangler, except that y appears as a parameter. 

Displacement surface. The displacement of the external flow streamlines by the 
retarding action of the boundary layer is an important effect. In two-dimensional flow 
the definition of displacement thickness, which is a measure of this effect, is straight- 
forward and there are several equivalent definitions. In three-dimensional flow, if one 
proceeds in strict analogy to the two-dimensional case, two displacement thicknesses 


can be defined [9]. These are 6, and 6, , where 
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in which (u, v) and (U, V) are the physical components of the velocity in the z’ and x” 
directions in the boundary layer and the external flow respectively and h is “some 
location well outside the boundary layer.” (The coordinates x«* are again general co- 
ordinates, not the streamline coordinates of the previous section.) Moore refers to 
these lengths as characterizing mass flow defects and then, using these,* defines a dis- 
placement surface as follows. In a geodesic normal coordinate system, the displacement 
surface x* = A(2', x’) is an “impermeable surface which would deflect a nonviscous 
fluid in such a way as to produce a normal velocity (W) satisfying”’ 


W = w(r', 2°) at 2° = A(z', 2’), 


where A has the same meaning as above, and W is the external flow normal velocity. 
After some approximations Moore gives a differential equation for A [9]. 

Here, using a slight modification of one of the two-dimensional definitions, an equa- 
tion for A is obtained. This approach is quite different from Moore’s but, making just 
the boundary layer approximations, this equation reduces to that of Moore. Analysis 
of the original equation shows an interesting difference between the three-and two- 
dimensional cases. 

In the following derivation geodesic normal coordinates are used and all vectors are 


in physical components. It is convenient to use the suffix notation for the range 1, 2 
but we set z = 2°. First we dispose of some geometrical preliminaries. It can be shown 
that (see [5], for example) the metric tensor components gg, are quadratic functions of z 
Jap = Gag + 2bape + Cape’, (11) 
where a,, has already been defined and 
bas = (Ogas ‘0z)o/2, 
Can = (8°gap/d2")o/2. 
Also the following relation holds 
Cap = G- On50q, - (12) 


Specializing to orthogonal surface coordinates, by means of (12), Eqs. (11) can be written 


as a perfect square. In the notation of (5) 


H, =h. + lz, (13) 

W here 
l. = bea/Ne (no summation). (14) 
If on a surface z = c = constant we draw a simply closed curve, the tangent vector, A, 


and the normal vector, n, have the components 
\: (HA, dx'/ds, H, dx’/ds), nr 
(15) 
n: (—H, dx’/ds, H, dx'/ds), 
where s is the arc length along the curve. 


To define a displacement surface we consider the flux of fluid through a developable 
surface S formed by the normals to the surface z = 0, passing through a simply closed 
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Fic. 2 
curved K, on z = 0, see Fig. 2. This flux is computed for that part of S between z = 0 
and z = h. For a velocity distribution Q this flux is given by 
ah . 
Fo = / Q-ndS = | $ (Q.H, dx’ — Q,H, dz*) dz, (16) 
S 0 + K 
where K, is the trace of the developable surface on the surface z = c and 0 < ¢ < h. 
The displacement surface 
z= A(z’, x’) 


is a surface such that the flux through S, for 0 < z < h, is the same for the following two 


velocity distributions 


I. Q,=Q=0 for 0<z<A 


% = for A<z<h 
G,= V 
A. iil for O<2<h 
(Q. v. 
From (16) we set 
F, Py 


and this is the condition imposed to find A. Interchanging the order of integration and 
making use of (13) and (14), after which the line integrals are calculated for A, , this 
requirement yields 


. 


p (Gho dx Gh, dx ) 0, (17) 
Ke 
where 


G, UA[1l + (1,4/2h.)] + [ (u — U){1l + (zl,/h.)] dz, 
si (18) 


h 
ah 


G. = VA[1 + (1,4/2h,)] + | (vy — V)[1 + (l,/h,)] dz. 


70 
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The integral in (17) can be written as an integral over the area bounded by K, using 
the “‘surface divergence theorem” [12]. Since K, is arbitrary we then obtain the result 


div G = 0 (19) 
for the vector G: (G, , G,) where the operator div is the “surface divergence’’ [12], i.e. 
div G = [a(h.G,)/dx' + A(h,G,)/dx’|/hyhz 
or in the subtensor notation (19) can be written 


G,, = 0. 


a 


lor given velocity distributions I and II and a value of h, (19) is a differential equation 
for A. However, in boundary layer theory, h cannot be given a definite value. In fact, 
in accordance with the requirements of this theory, we must let h — ©. (Extrapolating 
from the two-dimensional case, the integrals should converge since u — U, v — V 
exponentially as x — ©.) 

Thus far we have made no approximations. We take (U, V) to be the (physical) 
velocity components, in the z' and 2° directions, of the external flow evaluated on the 
surface z = O and (u, v) the corresponding velocity components throughout the boundary 
layer. Applying a transformation of the type (3) to z and A and keeping only the lowest 
order terms in (18) gives the boundary layer approximation to G, and G, 


G, = UA + [ (u — U) dz = U(A — 4), 
sis (20) 


i Pam [ =~ Vid= Vl ~ @), 


using the definitions (10). With these expressions (20) for the vector G, (19) is ess-ntially 
Moore’s equation [9] for the displacement surface. Moore gives some examples of the 
computation of A for special kinds of flow. Here we give a simple example by specializing 
only the coordinate system. For the streamline coordinates of the previous section 
V 0, therefore G, Jo vudz and (19) becomes 


afh,U(A — 6,)]/ax’ = —A(h,G.)/dx” = F(z", 2’). 


Thus 


h,U(A — 6,) - | F(x’, x’) dx' + f(z’), 
where f(2°) is a constant of the integration with respect to x' and must be evaluated 
from the conditions of the flow. Note that in the streamline coordinates the product 
V6. is well defined but 6, alone is undefined. 

The method used above to obtain A is a modification of the corresponding two- 
dimensional definition of displacement thickness. In two dimensions the developable 
surface S is taken to be a plane perpendicular to the plane of flow and as a result no 
differential equation need be solved for the displacement thickness. The A determined 
from (19) using (20) differs from the 6 of (10) by a constant of integration which, in 
most cases of interest, is zero. This is discussed by Moore [9]. 

However, it is interesting to compare the general expressions (18) for the vector G 
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in the two- and three-dimensional cases. Let 2° = 0 be the plane of flow for a two- 
dimensional flow. Then v = V = 0 and it is easy to show that 1, = 0. Thus the terms 
that are neglected after making the boundary layer approximations to obtain (20) 
disappear automatically for the special case of two-dimensional flow. This would be 
important if one wanted to calculate approximations to a viscous flow beyond the 
classical boundary layer theory as Kuo [13] has done for two-dimensional flow past a 
flat plate. For two-dimensional flow the displacement thickness expression does not 
change for the higher order approximations whereas for three-dimensional flow, including 
would be 


axially symmetric flow, an expansion of the expressions (18) in powers of v 


necessary. 
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OPTIMUM vs. CORRELATION METHODS IN TRACKING 
RANDOM SIGNALS IN BACKGROUND NOISE* 
BY 
R. C. DAVIS** 


Hughes Tool Company, Aircraft Division, Culver City, California 


I. Summary. A method of target tracking used exclusively in some applications 
is that of space diversity reception. We consider a two-dimensional problem with two 
receivers, although in principle the techniques developed are applicable to a three- 
dimensional problem utilizing three or more receivers. With N receivers the compu- 
tational difficulties increase as N*. The tracking problem considered here is formulated 
as follows. During a time interval 0 < ¢ < T one observes random processes 2x,(¢) = 


s(t) + n,(t) and x, (t 


= s(t + r)) + n.(t), and one desires to estimate the unknown 
value of ry) . s(t), n,(4), and n.(f) are continuous, Gaussian, stationary processes with 
continuous, monotonoid covariance functions and with Es(t)n,(t!) = Es(t)nt’) = 0. 
n,(t) and n.(t) are stationarily cross-correlated with a cross-correlation function that is 
in general an asymmetric function of time delay. As a criterion for an optimum estimate, 


ri, , of t) , we use the one-—-commonly used in statistical literature, that 75 should have 


minimum variance about the true value 7, ; that is, E(r3 — 79)” be a minimum. It is 
shown that for small values of w,7) (where w, is the highest frequency in the tracking 
pass band) the normalized error variance E[(74 — t»)/to]’ is identical to the reciprocal 


of the output signal-to-noise ratio—the criterion of common use in engineering literature 
on correlation methods. Using the Cramér-Rao inequality, an explicit expression for 
the minimum variance attainable by any estimate whatsoever is established. Moreover 
we arrive at the interesting conclusion that even with the use of optimum pre-detection 
filtering in a correlation system, the variance of the correlator estimate of 7, is always 
greater than the minimum variance given by the Cramér- Rao inequality. Gains in output 
signal-to-noise ratio obtained by an optimum system over a correlator vary in accordance 
with the statistical properties of signal and noise backgrounds. In general the gains are 
greater with asymmetrical cross-correlation between n,(/) and n,(¢) than otherwise. 
From a practical viewpoint gains obtainable over correlation methods become important 
for large values of the time-bandwidth product. This is due to the fact that under this 
condition the maximum likelihood estimate of 7) has an error variance that is approxi- 
mately equal to the minimum variance obtainable by any estimate. For small values 
of the time-bandwidth product no theory exists about the variance of the maximum 
likelihood estimate, since no efficient estimate of 7, exists. For this case, a worthwhile 
experimental study would be to compare the variance of the maximum likelihood 
estimate with the minimum variance calculated from the Cramér-Rao inequality. 
Finally it is shown how the maximum likelihood estimate is constructed and noted 
that the construction involves a considerable number of operations performed simul- 
taneously during the time interval 7. In order to simplify the system, an approximate 
maximum likelihood estimate is obtained that is optimum only ‘on target”, and its 


properties are discussed. 


*Received March 6, 1956. 
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II. Introduction. Heretofore problems involving several Gaussian stationary 
processes (such as the one outlined above) have been approached by the use of correlation 
methods that have been used extensively in engineering literature. With these methods 
the goal is the natural one of maximizing the output signal-to-ndise ratio. On the other 
hand there is a voluminous statistical literature on methods of analyzing problems of 
estimation of unknown parameters of probability distributions. The purpose of this 
paper is to apply some of these methods to the problem outlined. Using the same criterion 
of optimum as the one used in correlation methods, our method is more powerful. This 
is to be expected, since the best statistical techniques make use of the joint probability 
distribution of the coordinates of the various random processes entering into the problem. 
On the other hand correlation methods make use only of one second moment of the 
joint probability distribution. Since gains are demonstrated for Gaussian distributions 
(using all second moments in a particular combination), the author surmises that in 
problems involving non-Gaussian processes very large gains over correlation methods 
are to be expected. Unfortunately these gains are not obtained without difficulties. The 
main reason that modern statistical techniques have not been used extensively in noise 
problems is the formidable difficulty in obtaining coordinate systems that have the 
property that any finite number of coordinates are stochastically independent. In the 
case of Gaussian processes the central difficulty is that of matrix inversion. In problems 
involving a single random process, an orthogonal decomposition is possible by a method 
described by Karhunen [1]. In the case of two random processes (uncorrelated or not) 
possessing different power spectrum shapes, there exists no single set of orthogonal 
eigenfunctions that simultaneously yield orthogonal decompositions of both random 
processes. However, as we show later, in calculating the lower bound given by the 
Cramér-Rao inequality, we are able to avoid this difficulty. In order to obtain an explicit 
expression for the maximum likelihood estimate, however, we are forced to use much 
stronger assumptions on the statistical characteristics of signal and noises. 

III. Derivation of the Cramér-Rao lower bound. Consider the following problem. 
During the time interval 0 < ¢ < T we observe the random processes x,(#) = s(t) + 
n,(t) and z,(t) = s(t + ro) + n2(é). To is a fixed but unknown parameter. The problem 
is to determine the minimum variance of any unbiased estimate of 7, . We make the 
following assumptions concerning the characteristics of s(t), n,(¢), and n,(t). 

Assumption 1. s(t), n,(t), and n2(¢) are continuous, stationary, Gaussian processes, 
each with zero mean value, finite variance, and each possessing a monotonoid and con- 


tinuous covariance function. 
Assumption 2. s(t) is uncorrelated with n,(¢) and also with n,(?). 
Assumption 3. n,(t) and n.(t) are stationarily cross-correlated with a monotonoid 
and continuous cross-covariance function. 


It follows from Assumption 3 and a result due to Cramér [2, 


pio(t) = En,(t)n.(t + 7), 


. 227] that if we write 


then 


Pit) = | [cos wr dd,.(w) + sin wr dy;2(w)]. (3.1) 


The spectral functions ¢,.(w) and y,2(w) are of bounded variation in (0, ~). Each func- 
tion is the sum of an absolutely continuous function and a discontinuous function 
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(corresponding to possible discrete line frequencies in the spectrum). If no lines exist 
in the cross spectrum, (3.1) reduces to 


Pix(t) = | [di2(w) Coswr + Yi2(w) sin wr] dw. (3.1a) 
Similarly we denote by ¢,(w), ¢,(w), and ¢2(w) the power spectra of s(t), n,(¢), and n2(t) 
respectively. Each power spectrum is related to its corresponding covariance function, 


p(r), by the Wiener-Khintchine relation 


p(t) = | cos wt dow), (3.2) 
Jo 
where ¢(w) is of bounded variation in (0, ©). The first step consists in finding coordinate 
systems for x,(¢) and x,(¢) so that we can obtain an expression for the probability of 
realizing a given combination of 2x,(¢) and z,(¢). Since any practical system possesses a 
finite time-bandwidth product, we shall see that we require only a finite number of 
coordinates. To obtain the coordinates we use a slight extension of a method given by 
Root and Pitcher [3, p. 314]. It follows from Assumption 1 that we can expand p,(r) 
in a Fourier series valid in the interval —7 < +r < T [see 4, Chap. 6]. This gives 


p,(7) = . 6x COS woKr, 
Ke=l 


where w) = 2/7. Then there exists a process S(¢) defined by 
S(t) = Lim. } fax cos woKt + bx sin wK Et], (3.3) 
n-+@ K=1 


where the ax , bx are Gaussian variables satisfying 
Ea, = Ebxk = 0, for K = 1,2,--: 
Eajag = Eb;be = 8)x9« , Ea;bx = 0, j=1,2,---, K =1,2,-:-, 


and the process S(¢) has the same multi-variate distributions as s(¢). (This follows from 
the easily verified fact that ES()S(t + 7) = Es(t)s(t + 7) for every r < T.) Similarly, 
for the processes n,(¢) and n,(¢) we write 


p,(r) = ; ax COS w Kr, 
K=1 

p(T) = p> Bx COS w Kr, 
K-1 


pit) = >> [ye cos wKr + ex sin wKr], 


os 


> 


for —T < r < T, and correspondingly there exist processes N,(¢) and N,(¢) possessing 
the same joint multi-variate distributions as n,(¢) and n.(t). We define 


N,(t) = Lim. >> [ux cos wKt + vx sin w Kd), (3.4a) 
now K=1 

N,(t) = 1.i.m. Ye [px COS woKt + gx sin w Ki}, (3.4b) 
n-—-@ K-11 
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for 0 < i < T. Then we have 
Euz = Evi = ax, Ep, = Eax = Bx 
Euxpr = Evedx = Yx ; Eurdx = —Evepr = €x 
tor K = 1, 2, 
Moreover, every amplitude at frequency jzr/T is uncorrelated with every amplitude 
: | a I 
at frequency Kx2/T for 7 # K. 
From the above it follows that there exist Gaussian processes Y,(¢) and X,(/) which 
possess the same joint multi-variate distributions as x,(¢) and 2x,(¢). Clearly 
X,(t) = S()b + N,( = |.i.m. Zz [Ax CoS wKt + Bx sin w Kt], (3.5a) 
n-+ K=1 


Xe S(t + ro) + N.(t) = L.i.m. t B [Cx coswoKt + Dg sin w Kt), (3 .5b) 


for 0 < t< T. We have 


A, = Qn + Ux, EAx = 0x tax, 
By = bk + vx, EB = 0x + ax 
Ce = Ax COSwoKto + De SIN woK Tt, + Px ; ECxz = Ox + Bx 
D, —Arx sin woK tT, + be COS woK tT) + dk » EDs = 6.+ Br. 


The mixed covariances are as follows: 

EFAxBx = ECxDx = 0, 

BAC r Ox COS wk to + x ; 

EAxDx — Ox SiN wK > + €x ; 

EB .Cr = 0x SiN wK Tt, — €x , 

EB, Dre = Ox COSwK to + ¥x . 
Finally every amplitude at frequency 7j/T is uncorrelated with every amplitude at 
frequency *K/T, 7 # K. We are now ready to obtain the likelihood function of the 
first 4n amplitudes Ax , Be, Cx, De , for K = 1, 2, --- , n. Clearly this is the product 


of the likelihood function for each quadruple (Ax , Be , Cx , Dx). First we invert the 
4 & 4 moment matrix of (Ax , Be , Ce , De). The moment matrix is 


a 0 Ox COS WoK tT + Te —O9« SiNa KT) + €x 
0 Ox + ak Ox SIN woK 7, Ex Ox COS w KT, + YK 
Ox COS WK, + YK Ox SiINwaKty — «€ i. + £ 0 
—6x, sinwKry + €x Ox COSaKT, + 7 0 6. + B, 


The moment matrix can be inverted by the direct method of calculating the sixteen 
cofactors. A much simpler method is to partition the matrix into four 2 X 2 matrices 
and invert [5, p. 112]. The determinant of the matrix, denoted by Xx , is given by 


xe = (0- + ax)(Oc + Bre) — (Ox COS aK) + Yr) — (— 6x sin wKt, + €x). (3.6) 
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The inverse of the moment matrix becomes 








Ox + By 0 — 6x COS@Kt. — Ye Ox Sin wK tr) — éx| 
Nx Xx Ax 
Ox a Br = Ox sin woK To of — = Ox cos woK To -— oe 
0 eT Pe —<« x le 
XK AK AK 
0x cosanK'r, — 7x —OxSinwoKre+ ec Ox tag , 
XK AK AK 
Ox SINw KT) — Ex —Ox COS@K TT) — Yx Oc + ag 
Vx CIO Goes to F 0 Tx _T OE 
rz AK Xx 
Denoting by fx(Ax , Bx , Cx , Dx) the joint probability density of (Ax , Bx , Cx , Dx), 


we see [6, p. 311] that if we denote moreover by g, the joint probability density of the 


tn Gaussian variables Ax , Be , Cx , De with K = 1, 2, --- , n, that 


= f(A, ,B,,C,, Dif Az, Bs, C2, De) -** f(An ,B. , 


, 


ee Se 


: , IlIel,, . 
g. = (2m) TT dx! exp — 5 DES (0x + Bx) An + Br) 


+ (0¢ + axc(Ce + Dk) — 2(0x Cos wpK to + yx) AcCx (3.7) 


2 m 6x sin woK To -f- ex) AxDx ot 2(— Ox sin woK To 4 eB Cr 
— 2 Ox COs wo K To + Vx) Br Dx} ° 


ure now able to compute the lower bound for all unbiased estimates of 7, . In the 


We 
case of the multi-variate Gaussian distribution the necessary regularity conditions 
(loc. cit. p. 479) are satisfied, so that we can write 


OT»o 


E(r5 — 1) > |B 2. In a) | ‘ (3.8) 


For the problem considered here, we establish the fact that 


B| In 9.) = -(S, In tn), 


OT» 


and since the latter expected value is easier to calculate, we use the form 


- ,o Ing, | 
E(t — to) 2 | -2 no: | . (3.8a) 
OTs 


is established in the following manner. Consider the identity 
{oa J1 dq 8 In g,\ 
BA ae In ) = h\— ) _ z(* ——=)|. 
OT) g \9. OT, OTo 


Consider also the identity 


Equation (3.8a 


/ 


z(* 292) = - o p( 238-22), 
Jn OT OT) OTo 
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Hence (3.8a) is established if we show that (0 In g,/07)) vanishes identically in 7, . The 
reader can convince himself that (0 In g,/d7r,)) vanishes identically. Hence (3.8a) is 
established. 


From (3.7) we have 


3 eT | = f . 0° é,+ 8B 
S ) i ‘ 
3 In In 972 In X > > (A; t+ #3) ( ) 


OT 1 OF 
1 (C? + D?) 2; (= t #:) “442 (#: TO eette 3 ") 
OT) A; OTs Ay 


OTs dA, 
L ORK 0° ( 8; SIN wo; To + € 
2B ie 
OTe A 


2B.D, 2,(* pert 3 vw) 
T 


Taking the expected value of each side of the above equation and performing a tedious 
reduction, we arrive finally at a fairly simple expression. To emphasize the simplicity 


we define the following quantities: 


p 8. COS Wott. 4 Ves qi — 6, sin Wot To +é;, 
" Op ’ Od, 
p al qi = =. 
OT OTo 


In terms of these quantities, (3.8a) becomes 


ae gs wg os Ba a re 
iy? 3 < 25 \ (pi + gi) ry (PP. { 0a) (3.8b) 


a form involving the basic covariances and frequencies. 


Finally we reduce this t 


We have 
Pi + Qi = wot 8; , 
pp. + 4:9 Wt, (Y; SIN Wott, + €; COS WotTo). 


From (3.6) we obtain 


A O(a; + B 27; COS WotTy + Ze; SIN Wott.) + a8; — Yi — E& 


Denoting by w, the frequency woz, we obtain 
| , yY 6: 
Ne — -, “—~ 0 (a, + B 27 COS W;T) + Ze; SIN w;T)) + a,8; — Vi — 6 ies 
(3.8¢) 


w,0,(y; SIN Ww; To + €; COS w;To) 


A; 
}? 


—_ 1 0(a + B; — 2y,; COSw;To + Ze; SIN w;T)) + a8; — ¥; — €;) 


A glance at (3.8c) shows that the value of the right-hand side is independent of any 
linear pre-detection filtering performed common to both channels. Hence the theoretical 
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lower bound is independent of the bandwidth of the receivers in the two channels if 
they have identical amplitude-frequency responses. In this case (3.8c) becomes 


I ee — ; - w, 6° . 
K( 7, r,)* ~ = O(a; + B; — 2y; COS w;To + Ze; SIN w;T>) + a8; — v, —€, (3.92) 


. w,9,(y; SIN wT) + €; COS W; To) 


1 > 5——3? . 
— \0.(a; + B; — 2y; COS w;To + Ze; SIN w;To) + a;B; — ¥i — és 


Although the theoretical lower bound is given by (3.9a), we shall see later that this 
bound is never attainable by any practical system because of limitations of bandwidth. 
However, it will be possible to approximate the lower bound given by summing from 
i V, to. V, , where VN, — N, = [27'W], W being the receiver bandwidth, and 
|x] denoting the largest positive integer less than or equal to x. Finally we note that 
divergence of the right-hand side above implies that the lower bound is zero and hence 
that any estimate attaining this bound converges in mean square to the true value 75 . 
Later we shall see that for large values of the time-bandwidth product 7'W, the variance 
of the maximum likelihood estimate approximates the lower bound. 

IV. Performance of the finite time correlator. Although there is an extensive 
literature on correlation methods, there appears to be no discussion in the literature— 
at least for random signals——of the correlator that maximizes the output signal-to-noise 
ratio. Hence we give a brief exposition of this here. In order to facilitate a comparison 
of output signal-to-noise ratio with the lower bound of (3.9), we use the expansions 
(3.5a) and (3.5b). If we did not wish to make the comparison, it would be more natural 
to expand in terms of radian frequencies 2xn/T instead of the set mn/T that was required 
to obtain the lower bound. Since the Gaussian processes X,(/) and X,(¢) possess the 
same multi-variate distributions as x,(¢) and x,(t), we will obtain the correct value for 
the maximum output signal-to-noise ratio. The following is a block diagram of a cor- 


relator: 


x(t) | RECEIVER | 90° PHASE SHAPING 
AMPLIFIER SHIFTER | FILTER 
Mal Low-pass |¢o!t) 
__ FILTER 





x2(t) | RECEIVER | SHAPING | 
AMPLIFIER FILTER : 


In some correlators the 90° phase difference is obtained by using a differentiator in 
one channel, but this merely contributes a factor to the amplitude-frequency response 
of the shaping filter. To summarize, the phase of the signal in one channel is shifted 
90°—at each frequency —with respect to the other signal, each signal is filtered in the 
same manner by shaping filters, then multiplied together and low-passed to yield the 
output e)(¢). We shall determine the optimum pre-detection and post-detection filtering 
in order to maximize the output signal-to-noise ratio. Since we observe X,(¢) and X,(é) 
only during the time interval 0 < ¢ < 7, it will turn out that the optimum form of post- 
detection filtering is finite time averaging of the multiplier output. Hence we merely 
have to determine the proper form of pre-detection filtering. 

If we denote by Y,(¢) the process obtained by shifting the phase of each frequency 
component of X,(/) by 90°, we consider the processes 


r 
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X(t) = D> (Ax cos wxt + Bx sin wxt) (4.1a) 
K=1 
Y(t) = >> (—Cxsin wxt + Dx cos wxt) (4.1b) 
K=1 


for 0 < ¢ < T. Since we are interested only in the amplitude vs. frequency response of 
the shaping filter, for convenience we include the effect of this filtering in the Ax , Bx , 


Cx, and Dx . Denoting the amplitude response at frequency wx by fx’, we have 
Ax = (fx)'"(ax + ux), Bu = (fx)'"(bx + vx), 
Cre = (fx)'"(dx COS wet) + bx SiN WET, + Px), 
Dx = (fx)'"(—ax Sin wet + De COS wKTo + Qk); 


where the ax , bx , Ux , Ux , Dx , Vx are the random variables defined before in describing 


the processes X,(¢) and X,(¢). Clearly we have 


T co 


e= | XAOYA) dt = Od AC, + BBC, + AD; + €,B.D)), 
« i=l j 


a1) 
1 


a;; = 0, i+ Jj even, 
— + j odd, 
o; —w 
3, = 0 J 
—T/2, ;™= 3% 
y¥:;; = 0 all) 
T/2, $= 
«i = 0, i + j even 
2w 


——- i + j odd. 


Hence we see that 


Ee, = - dL E(-Be > + A,D,), 


= —T 2. f (0; sin w;7) — €;). 
1 1 
To obtain the output noise we first determine Ee, — (Ee,)°. We use the following identity 
valid for Gaussian variables each having mean value zero: 
E(X,X.X3X,) = EX,X,-EX,X, + EX,X;-EX.X, + EX,X,-EX.X, . 


Making use of the fact that any two amplitudes with different subscripts are uncorre- 


lated and using the values of a;; , 8;; , y;; , and ¢;; , we obtain after a tedious calculation 
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oO} 


en a ee 
i, (He) = > > fi) 26; SIN w@;To 
i eal | 


2) 
if 


+ O(a; + B; — 2, sin wT. — 2y; COS w;T) + a8; + € — ¥ 


she >» > ee { (ws + w°)(0, 4 «,)(0, + B,) (4.3) 


odd 
- Qa w,(0; COS w; To + ¥;)(0; COS w; To + Y¥;) 
2w;(8; sin w;T) — €,)(0; SiN w;t) — €;)}. 
Considering equation (4.2) we see that when p,.(—7T) # pi2(7), or equivalently that 


e; ~ 0 for at least one 7, that a portion of the d-c power is not useful signal. From (4.2) 


we write 
(Hey) | > £0; sin wire) + b> jes) — YY fi f;(O:€; sin wT) + 0;€; SiN w; To). 
I i=1 i=1 j=1 


Only the first term above represents useful signal power. Hence to determine noise 
power we must add the remaining terms. We denote by P,o/P,0 the output signal-to- 


noise ratio. We have 


P, E >> 7,0, sin wire] (4.4a) 


r’. = > £:{267 sin’ w; 7. 
= 1 1 
+ O(a; + B; — 2e; sin w;7) — 2y; COS w;T) + a8; + & — Yi} 
; ‘ (4.4b) 
+ T° > fue.) a p } FA (0;€; Sin wT) + 8;€; SIN w; To) 
i=] i=] j=1 
oi >» - Affi — ern 
(A@-—w)* 
odd 
The optimum correlator is determined by choosing the {f;}, 7 = 1, 2, --- ; ie., the 


amplitude frequency response of the pre-detection filter, so as to maximize P,o/P,0 , the 
output signal-to-noise ratio. Theoretically this is accomplished in the following manner. 
Since P,,» is clearly a positive definite quadratic form in the {f, , f;}, there exist orthogonal 
transformations [7, p. 13] that convert P,. to diagonal form; namely, > ee B°f; . Since 
P,, is the square of a linear form, it is converted by a linear transformation to the square 
of a linear form; namely, + Bel a,f,;)”. Hence 


ae (x xf.) | > ar] 


‘ i=l 
Applying Schwarz’ inequality to the above expression, we write 
=<} 2 -1 @ 
Pa/Pa = ( p a.fs) Pb af] < Diai/si. (4.5a) 
i=l i=l t=1 
Equality is attained in (4.5a) if and only if 
fi = Ka;/B; . (4.5b) 
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Hence (4.5b) vields the form of pre-detection filtering that maximizes the output signal- 
to-noise ratio. However, since f; > 0 for all 7, we see that the filter is physically realizable 
if and only if a; > 0 or a; < 0 for all ¢. If this is not the case, the maximum value of 
P,,./Py,. attainable is less than the value > Be , a,/B; . 

V. Comparison of optimum and correlation methods. It appears to be diffieult 
to make a general comparison of the extent to which the best correlator fails in having 
the optimum performance given by (3.9a). However, we can obtain some insight into 
the problem by considering a special case; namely, the situation arising when the noise 
cross-covariance function Is symmetric; 1.¢., py.(—7) pi.(r). First we show that for 
the practical situation in which ro is small that 75/1 (76 r,)” is identical with the 
output signal-to-noise ratio. For small values of wr, (where w, is the highest frequency 


in the receiver pass band), we see from (4.4a) that 


Qryv2 ae r2 2 
Po ~ ATA $00.) = Ker. 
1 
Moreover, ¢ Kr’, where 7’ is an unbiased estimate of 7, . Then 


’ ’ 2 r2yy_ 92 72 2 ’ , - 
P. Key — (Keo) K’Er K°7 KE (1 To) 
Hence, 
r To r’ ae i : 
a ro ; 2 ; o.1) 
F; NG To) 
Thus for small values of a7) , the output signal-to-noise ratio is approximately the 
reciprocal of the normalized error variance of 7’. 
Considering now (4.5a), we see that 
2 ryt 2 
a; ~ T'7 Ow; 
If now we approximate to Po by considering only the terms proportional to 7°, we have 


7 
pe = 120; sin’ w,;T) + O(a; + B 2y; COS w;T)) + a;B,; ete 


(This value of 8; will yield an output signal-to-noise ratio greater than the one actually 
attained, since the actual P,, is greater than the approximate one assumed.) Hence we 


obtain 


) 3 2 
Pro 25 Ou - (5.2) 
P — 20: sin’ w,T>) + 8O,(a r Dp 2y; COS W; T; t ap; Y, 

With the exception of the term 20; sin” w,7») in the denominator of (5.2), the above 


expression for P..o/P,. coincides with the first term of (3.9). In many cases of interest 
the first term of (8.9) is the dominant term; for example, for uncorrelated noises 
Pi2(7) 0 for all z—the second term in (3.9) vanishes. In those cases and to the extent 
that the approximation in (5.2) is valid, we can conclude that the output signal-to-noise 
ratio of the best correlator using finite time averaging as a post-detection filter is for 
all paractical purposes as great as any system whatsoever 

As an example of another type of correlator, consider one which uses an ideal low-pass 
filter with cutoff frequency 1/7 for post-detection, combined of course with an optimum 
form of pre-detectian filtering. It can be shown by using techniques similar to those used 


in the study of the finite time correlator that under the same assumptions that yielded 
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(5.2), symmetric noise cross-correlation and large values of the time-bandwidth product, 


the output signal-to-noise ratio is given approximately by 


aa 5 . w; 9; m ‘ (5.3) 
— 0,(a; + B; 27; COS W;T)) + a,B; — Yi 

Ilenee under these conditions the correlator using finite time averaging is 3 db 
superior to the correlator using an ideal low-pass filter with cutoff frequency 1/7’. 

VI. Considerations of time versus bandwidth. The practical considerations of 
receiver bandwidth necessitate a close examination of (3.9a). Since we observe the 
processes Y,(¢) and X,(¢) over the interval 0 < 4 < 7’, we obtain contributions to the 
output signal-to-noise ratio from arbitrarily high frequencies. When we discuss a method 
of extimating 7, so as to obtain optimum results, we shall see that the actual implemen- 
tation requires a receiver with infinite bandwidth and hence is not of practical interest. 
Before we finally restrict. ourselves to the practical case of finite bandwidth W, we wish 
to note an interesting result. An examination of the right-hand side of (3.9a) shows 
that the infinite series either converges or diverges depending upon the relative rates 
with which the signal and various noise spectra approach zero as the frequency approaches 
infinity. Hlenee we conclude that finite time observation using infinite bandwidth need 
not result in an infinite output signal-to-noise ratio except under special conditons. 
On the other hand consider the situation that prevails when the bandwidth is finite 
and restricted to the range (Nia , Now), where N, — N, = [27'W]. We wish to obtain 
the limiting form of (3.8c) as T’ —> ©. Since the interval between adjacent radian 
frequencies is r/7', we see that by multiplying and dividing the right-hand side of (3.8¢) 
by 2/7’, the sum approaches an integral, and we have for large 7, 


I 
Ki r ) 
4 fl w [O(w) }* dw etree 6.1) 
rT. O(w) [a(w) + B(w) — 2y(w) cos wr) + 2e(w) sin wro]-+-a(w)B(w) — [y()]’—[ew)]° >" 
17’ [ wO(w) [y(w) sin wro+e(w) COS wTo| ya 
: Kt 57 dw. 
s | A(w) la(w) + B(w) — 2y(w) COs wry + Ze(w) SiN wT] + a(w)B(w) — [y(w) |” — [e(w) | 


Hence we arrive at the following conclusion—with the obvious restriction that the 
tracking pass band includes a positive amount of signal spectral energy. With finite 
observation time and infinite bandwidth the output signal-to-noise ratio may or may 
not be infinite in accordance with the statistical character of signal and noise. On the 
other hand the combination of finite bandwidth and infinite observation time always 
results in an infinite output signal-to-noise ratio. 

VII. Asymptotic properties of the maximum likelihood estimate of r, for large 
values of time-bandwidth product. We emphasize the fact that in this section we 
restrict the analysis to the practical case of finite bandwidth. In order to accomplish 


this, we rewrite the lower bound for the variance as 


| iti \ w 0? 
E(x 7 ~ we Oia. + 5, Qy; COS w,;T)o + Ze; SIN w;T)) + a,B; — v: —€ (7.1) 
- 


a w,O,(y; SIN w; To + €; COS &; To) \ 
SX, (O(a; + By — 2y; cos wT) + Ze; SiN w;T)) + a,8B; — Yi — é , 
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V, and N, correspond to the lower and upper radian frequencies N,r/T and N.7/T, 
where as before VN, — N, = [27 W], W being the receiver bandwidth in cycles per second. 

We consider now an estimate of 7) which for large values of 27W has a variance 
that approaches the right-hand side of (7.1). This is the maximum likelihood estimate 
introduced by R. A. Fisher (for a discussion see Cramér’s book, loc. cit. pp. 498-504). 
Since the properties of the maximum likelihood estimate follow directly from properties 
of the likelihood function for the random processes x,(¢) and x,(¢), we use in lieu of this 
the likelihood function g, of Eq. (3.7), corresponding to the processes X,(¢) and X,(¢) 
defined by (3.5a) and (3.5b). As stated earlier, this is justified since the processes x, (/) 
and x,(¢) possess the same multivariate distributions as X,(¢) and X,(¢). It follows that 
the maximum likelihood estimate of 7) determined from the processes X,(¢) and X,(t) 
will have the same statistical properties as the corresponding estimate determined from 
x,(t) and x,(t). We emphasize now, however, that since X,(¢) and X,(¢) are not the 
processes observed, one is still left with the problem of constructing the estimate from 
the observed processes x,(¢) and x,(¢). This question will be considered later. 

There is very little material in statistical literature on properties of the maximum 
likelihood estimate of an unknown parameter in a continuous random process. Grenander 
has shown [8, p. 255-257] that for a process of generalized Markoff type, the maximum 
likelihood equation has a root that is consistent and asymptotically efficient. More 
recently P. Whittle [9] has shown that under certain assumptions on the spectrum of a 
single Gaussian process, the maximum likelihood equation has a root that possesses the 
usual optimum properties. We are confronted, however, with estimating a parameter 
of the joint distribution of two processes. Wald [10] has investigated the asymptotic 
properties of the maximum likelihood estimate of an unknown parameter of a discrete 
stochastic process. We are able to use a slight extension of his method to obtain our 
results. The maximum likelihood equation is given by 


ar, In gy = O, (7.2 


where 


N = [2TW]. 


We add another assumption to Assumptions 1-3. 

Assumption 4. The power spectral density of the signal, s(/), has only a finite number 
of zeros in the pass band N,2/T to N.r/T. This assumption is not any practical re- 
striction on the generality of the analysis, and is as a matter of fact stronger than 
required. We denote by Cy(7.) the reciprocal of the right-hand side of (7.1). A sequence 
{rv}, (N = 1, 2, --- , ad inf.) of estimates of 7, is said to be asymptotically efficient in 
the wide sense if the mean of [C'y(1,)]'”? (7% — 70) is zero and the variance of [Cy(r)]'? 
(rx — 7.) is 1 in the limit as N — o. Instead of the joint probability density 


gy(A, , By ’ C; 9 D, ; As 9 B, ‘ C. 9 Da a Ss Aa , By ) Cn ’ Dy) 


which we have, Wald considers a probability density py(X, , X2 , , Xy , To) defined 
for all 7) in some non-degenerate interval A on the real axis. The following four con- 


ditions are assumed to hold: 
Condition 1. The derivatives d‘py/dr, (i = 1, 2, 3) exist for all 7) in A and for all 


samples (x, , %2, ‘+: , Zy) except perhaps for a set of measure zero. We have furthermore, 
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° x | af | 
fos flab. | 29%) dey dey +++ dey << = = 1,2). 


rotA | O75 | 





Condition 2. For any 7» in A, limy... Cy(to.) = ©. 

Condition 3. For any 7» in A the standard deviation of 0° log p,/dr¢ divided by 
the expected value of d’ logp,/d75 (both computed uder the assumption that 7» is true) 
converges to zero as N > o, 

Condition 4. There exists a positive 6 such that for any 7, in A the expression 


[| 
is a bounded function of N for 74 in the interval | 73 — 7>o| < 6. 
Conditions 1 and 4 are clearly satisfied, since gy(A, , B,,C,,D,,-:-,Aw, By, Cw , Dw) 


EB d° In p,(tr y Ze y *** » Ew y 70) 
Culte) 


l.u.b. | = 
43 
re? OTs 





is of the form 


Ky(r ) exp {—Qy(A, , B, ’ C; ’ dD, poe Ay , By ’ Cy ’ Dy/t)}; 


VJ 


where Qy is a positive definite quadratic form in the variables Ax , Bx , Cx , Dx , K = 
N,, «++, N.. Condition 2 was shown to be satisfied in Sec. VI, since the output signal- 
to-noise ratio approaches infinity as 7’ — ©. Condition 3 can be shown to be satisfied 


in the following manner. First we have 


42 i=N, 92 
- ae = 
“Ingyv = > san fA; , B; , C; , Dd), (7.3) 
OT) : ijn, OTo . 
where as usual V, — N, = N. Since this is a sum of N uncorrelated random variables, 
we have 
2 a. i=N 2 2 4? 2 
a , Oo” 4 ; 0 > i 0 if , 2 
P( “In gy) - (B ; In ) = oS, In ) _ (e —; In J} <N 
\OT5 ‘9g ) OTo gn >» OTo J OTo J ania 


where 
Moreover we have 


W here 


0 


(oa Sade 
6 = inf Be 3 In fi). 
; 0 


From Assumption 4, 6 > 0. Hence, 

( 42 \ 2 42 2) 1/2 2 

ee Oo 2 0 : 0 -lar-1/2 , 

¢ B(: 3 In gv) (EZ —; In a) K( 5 In aw) <o6 JN —-o as Noo, 
OT) P \ OT) OT 

Since Wald’s Conditions 1-4 are satisfied, his analysis is applicable directly to our 
problem. Then we arrive at the result that Eq. (7.2) has a root that is a consistent estimate 
of r. . Furthermore any root of (7.2) that is a consistent estimate of 7. is also asymp- 


totically efficient in the wide sense. 
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We have determined the asymptotic properties of the maximum likelihood estimate 
by analyzing the joint distribution of the random processes X,(¢) and X,(t). In order 
to determine the maximum likelihood estimate explicitly we require the joint distribu- 
tion of the observed processes x,(¢) and x.(/). The determination of this distribution 
under the weak Assumptions 1-3 requires the inversion of an N X N matrix, with NV = 
{27 W]. Since the matrix is nonsingular, this inversion can be accomplished in principle. 
In order to see the required filtering operations on 2,(¢) and x,(¢) more clearly, we make 
an additional assumption. 


Assumption 5. Each of the covariances p,(r), p,(r), and p2(r) is a periodic function 
p pi P2 


of 7 with period 7’. 
Clearly it follows from Assumption 3 that p,,(r) is also periodic with period 7. As 
Root and Pitcher have shown (loc. cit., p. 313), it follows then that s(t), n,(¢), and n.(t) 


possess Fourier expansions in (0, 7’) with pairwise orthogonal amplitudes. We write 


s(t) l.i.m. 2. [ax cos QKt + be cos QKtI, (7.4a) 
n-> K=1 
n,(t) lim. >> [ux cos 2Kt + vg cos %Ke], (7 .4b) 
n— K=1 
n 
n(t) = lism. Zz. [px cos 2,Kt + qx sin 2 Kit], (7.4¢) 
n+ K=1 
for < ¢ < 7. 
We note that 2, = 27/7’, whereas the processes S(t), N,(¢), and N.(t) were expanded 
into components with frequencies being multiples of w = 2/7. Then we write 
a(t) = s(t) + n,(t) = lim. DO [Ax cos 2,Kt + Bg sin Kt], 
n—- © K=1 


' 


x(t) = s(t + 7) + n.(t) = 1.i.m. 2. [Cx cos QKE + Dx sin Q Kt], 
n-+ K-=1 


for 0: < i < T. 

To save useless reiteration we use the same notation for the second moments of 
Ax, Be , Cx , Dx as we did for the corresponding amplitudes of the processes X,(¢) and 
X,(t). Denoting by h,(A, , B, , C, , D, , --- , An, B, , C, , D,) the joint probability 
density of the Ax , Bx , Cx , Dx , the maximum likelihood equation becomes 


0 —" 1 2 ea ( a », (0; + B;) 
a Sigg Inds — 5 ((A? + Bye 
OTs =o 2, OT. . 2 OT r» + & d; 
+ (C2 + D) PE _ 94,0, + B,D) A ein ty) 7s) 


— sin Q,t7, ) 
~ (A.D, — B.C, (— 6, sin oe + a = 


Let us consider what filtering operations are required to obtain the terms in the likelihood 
equation. The basic tool required is Parseval’s equation for Fourier series. This states 


that if we have the Fourier series 
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y(t) = ¥ (a, co Ss om 4 b, sin art), 


y2(t) -” PB (c COs “ fs d,s in 7"), 
tel 
forO0 <?¢ < T, then 
9 aT > 
- | y,(t)y(t) dt = > (ae, + b,d,). 
#0 t=1 
Examining the term 
i No 
>> (Ai + pa( ot bs), 


we see that if x,(/) is passed through a linear filter with gain [(@, + 8,;)/,]'” at frequency 
0,2, yielding an output y,(¢), then 


~T 


- (A? + B) (24 ca Bs) = =f [y(F dt. 


An analogous operation is performed on 2x,(¢) to obtain the term 


N \ 
i a ; 6, 
> C+ Di a ). 
‘<M, Aj 
lo obtain the term 
N . \ 
—* f i] 208 Q) 0 £ 
S (AL, + B.D)| cet ts), 
7, r; 
it is neeessary to pass 2,(4) and 2,(¢) through identical linear filters with gain 
[(@, cos Qito + y,)/\,]'~ obtaining outputs Z,(¢) and Z,(t). Then 


a -08 } 9/7 } 
> AC. + B.D,” CO! Seite + 9s) = T | Z,(t)Z,(t) at. 


In order to obtain the term 


_ 6, sin ito + & -t) 
24 (A;D, — B,C (= =" | 


i 


we note that we must first pass x,(¢) through a 90° phase shifter, obtaining at the output 
p £ I ’ £ ] 


t=N, 
z(t) = >» (—C; sin Qt + D; cos Qoil). 
Ni 


Then x, (¢) and x-(¢) are passed through identical linear filters with gain 


(= 0, sin Qt +)" 
A; ‘ 


producing processes €,(t) and e,(t). Then 


i=N, = a. Os 9 ~T 
> (AD; = B.C) = sin Qin» + 41) == | e(t)e,(t) dt. 
A; 7 /0 


iN, 
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It is interesting to note that the last term in the likelihood equation (7.5) requires 
exactly the same operations as the finite time correlator; namely, 90° phase shifting of 
one process, linear filtering, multiplication, and finite time averaging. 

Brief consideration of Eq. (7.5) shows that it cannot be solved explicitly for 7} in 
terms of the various filtered quantities. One possible approach is to evaluate In hy for a 
sequence of values of 7) and determine its absolute maximum. Clearly this involves a 
large number of operations in comparison to those required by a simple correlation 
system. In the next section we discuss a method for obtaining an approximate maximum 
likelihood estimate. 

VIII. An approximation to the maximum likelihood estimate. In order to avoid 
the obvious difficulties in solving (7.5), we discuss an approximation to the maximum 
likelihood estimate. The method is based upon the fact that in the process of tracking 


one is always striving to achieve the condition 7, = 0. Hence consider 0/079 In hy evaluated 


at rt, = 0. This condition prescribes the filtering operations for x,(¢) and x,(¢) for the 
“on target” situation. Let us denote by (0/07, In hy) the expected value of 0/079 In hy 
computed when 7 is true. We have shown in See. III that #,,(0/d75 In hy) = 0. Hence 


this estimate has no average d-c bias due to noise for the ‘“‘on target’? condition. More- 
over, the output signal-to-noise ratio for small deviations from 7, = 0 attains the maxi- 
mum value possible. As the true value of 7) continues to deviate from zero, the output 
signal-to-noise ratio deteriorates below the maximum attainable. 
It is of practical interest to investigate the “steering pattern” of this approximate 
.o . After a few algebraic manipulations, we obtain 
oer 


(a ‘eX? Dy .62 J, 7 - 
B= In hy) = p> 52 * 10 + €;) sin w;r + 2€,(0; + y,) sin” “of. 


estimate; i.e., H,.(0/dr_ In hy),, 


Hence we see that the average d-c output consists of the sum of a symmetric term and 
an anti-symmetric term. For background noise with symmetric cross-correlation (¢ 0 


for every 7), the output consists solely of a symmetric term. 
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THE APPLICATION OF LIMIT ANALYSIS TO THE DETERMINATION 
OF THE STRENGTH OF BUTT JOINTS* 


BY 
R. T. SHIELD 


Brown University 


Summary. ‘The technique of limit analysis is applied to determine upper and lower 
bounds for the tensile strength of a butt joint consisting of a thin layer of adhesive 
joining the parallel flats of two rigid adherends. The adhesive is assumed to be an elastic- 
perfectly plastic material which yields when the maximum shear stress reaches a critical 
value. The methods used apply to any joint with a convex area of cross-section. Par- 
ticular application is made to joints whose cross-sections are circular, rectangular, or a 
polygon circumscribed about a circle. 

1. Introduction. It is known experimentally (see [1]**, for example) that the 
strength of a butt joint of a hardened adhesive between two rigid adherends is inversely 
proportional to the thickness of the adhesive layer when the layer is thin. Various 
theories have been suggested to explain this phenomenon and critical reviews of these 
theories can be found elsewhere [2, 3]. The theory of the present paper rests on the 
assumption that the adhesive can be represented with sufficient accuracy by an elastic- 
plastic material which yields under constant stress when the maximum shear stress 
reaches a certain critical value. Under this simplifying assumption, the problem of the 
determination of the stress distribution in a butt joint under tension is still one of con- 
siderable difficulty. However the assumption enables the technique of limit analysis 
[4, 5, 6] to be used to predict the strength of the joint. In Secs. 3 and 4 below, particular 
application is made to joints whose sections are rectangular, circular, or a polygon 
circumscribed about a circle, but the analysis can be modified so that it applies to any 
joint with a convex area of cross-section. An approximate analysis for a circular joint 
has been made previously by Kachanov [7]. 

The related problem of the plastic compression of a layer between two rough plates 
has also been treated by approximate analysis [8, 9]. Since the yield condition used 
here for the adhesive is independent of hydrostatic pressure, reversal of the stresses in 
the analysis below gives results which apply to the related problem, and these results 
agree with the corresponding results in [8, 9]. 

2. Limit analysis and statement of the problem. For an assemblage of rigid and 
of elastic-perfectly plastic bodies, the term collapse will be used here to describe con- 
ditions under which plastic flow would occur under constant loads if the accompanying 
changes in the geometry of the assemblage were neglected. Assuming that the boundary 
conditions are of the stress type (i.e. each component, 7’, , 7, or T, , of the surface 
traction is given except where the corresponding velocity component, v, , v, or v, , or 
the corresponding relative velocity component at the interface of an assemblage is 
prescribed to be zero), the limit analysis theorems [4, 5, 6] provide a method for determin- 


*Received March 9, 1956. The results presented in this paper were obtained in the course of research 
sponsored by the Office of Naval Research under Contract Nonr-562(10) with Brown University. 
**Numbers in square brackets refer to the Bibliography at the end of the paper. 
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ing whether or not collapse will occur at any stage in a given loading program. The 
limit analysis theorems of Ref. 6 which will be used below are outlined briefly here. 

A stress field is said to be statically admissible if the stress field is in equilibrium 
and satisfies the stress boundary conditions. The stress fields employed in the following 
consist of regions of constant stress separated by surfaces across which the stresses are 
discontinuous. These fields will be in equilibrium if the expression 


OT2Ns TT TeyNy TH Tass 5 (1) 


and each of the two similar expressions, is continuous across a surface of stress dis- 


continuity, n, , n, and n, being the components of the unit normal to the surface. A 


stress field is said to be safe if the yield condition is nowhere violated. 
From any velocity field v, , v, , v, , strain rates can be derived in the usual way; 
Ov, Ov, Ov, 


 ~. Re oe _— »* . 


Oy Ox 
If these strain rates are treated as purely plastic strain rates, the internal rate of dissi- 
pation of energy can be calculated, provided that the strain rates satisfy any condition 
imposed by the yield condition. Here the yield condition is assumed to be independent 
of the hydrostatic component of the stress tensor and imposes the incompressibility 


condition, 


é+¢, +e, = 0, (3) 


on plastic strain rates. A velocity field is called kinematically admissible if it satisfies 
the incompressibility condition and satisfies the velocity boundary conditions. Such a 
velocity field is said to be a kinematically admissible state of collapse if the total internal 
rate of dissipation of energy is not greater than the rate at which the applied tractions 
do work on the velocities of their points of application. 

The following theorems have been formulated for the case of all surface tractions 
increasing in proportion: 

Theorem 1. Collapse will not occur until the largest values of the surface tractions 
are reached for which it is possible to find a safe statically admissible stress field. 
Theorem 2. Collapse will occur under the smallest values of the surface tractions for 
which it is possible to find a kinematically admissible state of collapse. 

Theorems 1 and 2 can be used to determine lower and upper bounds, respectively, 
for the collapse values of the surface tractions. 

For a material obeying Tresca’s yield criterion, plastic flow can occur under constant 
maximum shearing stress k, and safe statically admissible stress fields in the material 
must nowhere involve a shearing stress greater than k. Also, for this material, the internal 
rate of dissipation of energy per unit volume due to a plastic strain rate is 2k max | ¢€ |, 
where max | «| is the numerically largest principal component of the plastic strain rate 
[10]. Thus kinematically admissible collapse states must be such that 


/ (Tv, + Tp, + Tv) dS > | 2k max | «| aV, (4) 
Ss JV 


where S is the bounding surface and V the volume of the body or assemblage of bodies, 
and where max | ¢ | refers to the strain rates derived from the velocity field. Under the 
Tresca criterion, discontinuities in the tangential velocity across fixed surfaces are 
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permissible in kinematically admissible velocity fields. At a surface of discontinuity, 
energy is dissipated at the rate kAv per unit surface area, where Av is the magnitude of 
the relative change in velocity across the surface, and this must be taken into account 
in evaluating the total internal rate of dissipation of energy. 

For the purposes of this report, a butt joint is treated as a thin layer of adhesive 
joining the parallel flats of two rigid adherends. The free ends of the adherends are 
loaded by two equal and opposite tensile forces of magnitude 7’, having a common line 
of action which is normal to the joint and which passes through the centroid of the 
cross-sectional area of the joint. The adhesive is assumed to be an elastic-perfectly 
plastic material obeying Tresca’s yield condition. The tensile force T across the joint is 
increased until the collapse load is reached and the joint fails due to plastic yielding of 
the adhesive layer. The adherends and the adhesive layer constitute an assemblage of 
bodies to which the limit analysis theorems can be applied to determine the collapse 
value of the load 7. Statically admissible stress fields in the layer must be such that 
there is no shearing stress on the mid-plane of the layer because of symmetry, and also 
the resultant tensile force on the mid-plane must act through the centroid of the area 
of the mid-plane. At the interfaces of the adhesive and the adherends, the tractions can 
have any value compatible with safe statically admissible stress fields in the adhesive, 
since the adherends are rigid and since there is complete attachment between the adhesive 
and the adherends. Kinematically admissible velocity fields must be such that there is 
no relative motion between the adherends and adhesive immediately adjacent to the 
adherends, which move as rigid bodies. However, surfaces of tangential velocity dis- 
continuity between the adhesive and the adhesive immediately adjacent to the ad- 
herends are permissible. 

3. Lower bounds. In this section, a method is given for obtaining lower bounds 
for the collapse value of the tensile load T for butt joints which have a convex area of 
cross-section. The stress fields employed are derived from a stress field used previously 
to obtain a lower bound for the average indentation pressure in the plastic indentation 


of a layer by a circular punch [11]. 


! 

| (2+2/2)k J (2+/2)k 
a a —- =f F | A. 
i 
| 
Cc” 7 


’ 
——— 

















Fig. 1. Stress field for plane strain problem. 


Figure | shows the plane strain version of the stress fields used. In the figure, OA 
is the mid-plane of the adhesive layer of thickness 2h and BC is the upper surface of the 
lower adherend of width 2a. The field is symmetrical about the mid-plane 0A of the layer 
and about the central plane OC. Region ABD is unstressed and lines BD, DE, EF, ete. 
are lines of stress discontinuity. The tensile stress of amount 2k in region BDE parallel 
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to BD produces a vertical tension (2 + 2'”)k and a horizontal tension 2'”*k in region 
DEF, where k is the yield shear stress of the adhesive. The stress in region EFG is taken 
to be a hydrostatic tension of amount 2'’*k in order that a tension 2k parallel to EG 
can be superimposed in region EGH. This stress increases the stress across GJ to (2 + 
2°”*)k. By repeated application of this process the tension across the mid-plane is in- 
creased towards the center O. 

The extension of the field depicted in Fig. 1 to a butt joint of circular cross-section 
is shown in plan in Fig. 2. From strip elements of area on the conical surface through 
DE in Fig. 1, “legs’’ of material originate and carry tensions of amount 2k inclined at 
an angle of amount 2/4 to the horizontal. Thus DEE’D’ in Fig. 2 is a strip element of 














Fig. 2. Plan of stress field for circular joint. 


area on the conical surface through DE in Fig. 1 subtending an angle 60 at the axis OC. 
This strip element generates the “‘leg’’ DEBB’E’D’ which carries a tension 2k in the 
direction BD. The material vertically above the conical surface through DE is stressed 
by a vertical tension (2 + 2'*)k and equal radial and circumferential tensions 2'’*k. 
The material lying between the cylindrical surface through /F and the conical surface 
through EG is stressed by a hydrostatic tension 2'*k. The process of increasing the 
tension across the mid-plane towards the center of the punch is repeated, so that the 
annulus in the mid-plane bounded by the circles through G and J carries a normal 
stress (2 + 2°*)k and so on. Further, triangular elements of area in the mid-plane, 
such as B’D’K in the plan, do not lie vertically above “legs” of material carrying a 
tension 2k. The tension over these areas can be increased by adding a vertical tension 
2k in the triangular prisms of material below the areas. 

The field described nowhere violates the yield condition and is a statically admissible 
stress field, since it satisfies the symmetry conditions on the mid-plane. It therefore 
, where n 


provides a lower bound for the collapse load of the joint. When a/h = n2’ 
is an integer, the process of increasing the tension towards the center is carried out n 
times, and the average tension ¢, over the mid-plane is found to be given by 


t, is, 14 in  B\h 
1-2 °° + +{2°° — . (5) 
k 3h 3 
The product of ¢;, and the area ra’ of the joint is a lower bound for the collapse load 7’. 
The stress field outlined above can be modified so that it applies to any joint which 
has a convex area of cross-section. When the glued section is a rectangle with sides 2a, 
2b, where b > a, the average tension /, over the mid-plane is found to be given by 
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i, ue, ROR. Da = 2) h 
aa _ _— 4 
k ; zs sl 3h (3 2 *) - (2 af 6’ (6) 
for a/h = n2'’°, where as before 2h is the thickness of the adhesive layer and n is an 


integer. The value 4abt, is a lower bound for the collapse load T’. 

For a square cross-section, a = b and the value (6) reduces to the value (5). Indeed 
it can be shown that for any joint whose cross-section is a polygon circumscribed about 
a circle of radius a, the average tension ¢; over the joint provided by the stress field is 
given by (5) for a/h = n2"'””, 

An alternative stress field for a rectangular joint can be derived from a stress field 
used previously in [11] to obtain a lower bound for the average indentation pressure in 
the plastic indentation of a layer by a square punch. The field provides a value ¢; for 


the average tension over the joint given by 


t, la/3 1 2) h (2 a ‘) 1h? "i 
} 1+ 24(3_ ie vo 6b 2 + 3 ab? (4) 
for a/h | + 2'* + n2°*”, 2a and 2b (b > a) being the sides of the rectangle. The 


value (7) for ¢; is greater than the value (6) but the percentage difference between the 
values is small for large values of a/h. 

4. Upper bounds. In applying Theorem 2 of Sec. 2 to obtain upper bounds for the 
collapse value of the tensile load 7 across a butt joint, the procedure is as follows. A 
mode of deformation of the adhesive layer is assumed which is compatible with the 
incompressibility condition (8) and the condition that the adherends move as rigid 
bodies. The total rate of plastic dissipation of energy due to this velocity field is then 
calculated. With this kinematically admissible velocity field the rate at which work is 
done by tensile loads 7' applied to the adherends can be found, and the value, 7'y say, 
of T is chosen so that the rate at which the loads 7’, do work is equal to the rate of plastic 
dissipation of energy in the layer. Under the loads 7’, the velocity field is a kinematically 
admissible state of collapse and it follows from Theorem 2 that 7’, is an upper bound 
for the collapse value of the load 7’. 

In the velocity fields employed below, it is assumed that the adherends have equal 
and opposite velocities V normal to the joint so that the adherends separate at a velocity 
2V. The z-axis is taken normal to the joint area and the mid-plane of the adhesive layer 
is taken to be the plane z = 0. The velocity fields are symmetrical about the plane z = 0 
so that only the motion in the upper half of the layer, the region 0 < z < h, need be 
considered. The velocity component v, is zero on z = 0 because of symmetry and has 
the value V on the flat z = h of the upper adherend, since v, must be continuous across 

h. These conditions on v, are satisfied by assuming that 
vy, = V2/h, O<z2<h. (8) 
The velocity components v, and v, in the layer must then satisfy the incompressibility 
equation 
wal 4 Ov, =e if (9) 
Ox oy h 


and for simplicity v, and v, are assumed to be independent of z. With cylindrical polar 
co-ordinates (r, 6, z), Eq. (9) becomes 
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Ov, v, 1 Ov, V 10 
= ( ) 

or r r 00 i 


where v, and v, are the radial and circumferential velocity components, assumed to be 


functions of r and @ only. 
For the joint of circular cross-section we take 


v, — Vr/2h, Ve 0, oO < ¢ <a, 11) 


where the origin of the coordinate system is at the center of the middle section. On 


z = h, v, and v, are zero. Energy is dissipated in the layer due to the non-zero strain 


rates and also in the velocity discontinuity surface between the adhesive and the adhesive 


immediately adjacent to the flats z = + h of the adherends. A straightforward calcula- 


tion shows that the velocity field will be a kinematically admissible collapse state if the 





average tension over the mid-plane is greater than ty , where 


t, la 
eee 
k - t = a 


(12) 


a being the radius of the cross-section. By Theorem 2 the product of fy and the area 


ma of the cross-section is an upper bound for the collapse value of the load 7. The 
percentage difference between the upper bound (12) and the lower bound (5) for the 
average tension over the joint tends to zero as a/h tends to infinity. 

The simple velocity field (11) can also be used for non-circular cross-sections. For a 


joint with a square cross-section of side 2a, the field is a kinematically admissible collapse 


state if the average tension over the joint is not less than ¢, given by 


U = 2 + 0.383 2. 13) 
k h 
For very large values of a/h, the difference between the upper bound (13) and the lower 
bound (5) is 15 per cent of the lower bound. 
In order to improve the upper bound (13) for a square section, and also for applica- 
tion to other sections, more sophisticated velocity fields are required. With a view to 


constructing velocity fields for polygonal sections, consider a right prism of adhesive 


whose section is a right-angled triangle, triangle OAB in Fig. 3, and which is bounded 


by the planes z = 0, h. The coordinate system is taken as shown in the figure and 


v, is again given by Eq. (8). The components 7 
0 < z < h, the incompressibility condition (10) is satisfied by the expressions 


and v, are zero on z = h. In the region 
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Fic. 3. Velocity field in triangular region. 
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where f is a function of @ only. We shall assume that v, is zero on the sides 6 = 0 and 
6 = a of the prism since this will facilitate the fitting together of fields for a poly- 
gonal cross-section. The total rate D of dissipation of energy in the prism and in the 
discontinuity surface z = h can, in principle, be minimised with respect to the function 
f(@), subject to f(0) = 0, f(a) = a. This will not be done here, however. Taking f(6) 
to be given by a linear function of 6 gives the field (11). Assuming f(6) to be a quadratic 


function of @ gives 
6° r 
f(0) = Bo(1 _ *) += (15) 
a a 


and the associated rate D of dissipation of energy can be minimised with respect to the 
arbitrary constant B. The minimum value of D considered as a function of B is near the 
point B = 2 and taking B = 2 provides results sufficient for the purposes of this paper. 
Substituting (15) in (14) and setting B = 2 gives the components 


= =— vr (1 - 2), % = Ve i - $), (16) 
h a h a 


It should be noted that the value of v, on @ = 0 and 6 = a is independent of a. It is 
found that the total rate D of dissipation of energy is given by 


D fe 
= a 1+ M(a) + 3h K(a), (17) 
where 
bs 2 a a@ : 4 e : " 

K(a) = | sec’ 0 (1 ~ Ja + 0)” dé, (18) 
tan @ Jo a 

M(a) = ; [ sec’ 6 g(@) dé (19) 

: tana Jy . s 


1 a’ tan a@ is the cross-sectional area of the prism. In expression (19), 


¢ 2°]1/2 
g(0) = mast, | 4. + (1 - 2)" \ 
Qa Qa 


The values of K(a) and M(a@) for values of a between 0° and 90° are shown in the table 
to three decimal places. As a tends to zero through positive values, M(a) tan a tends 


and where A 
g(@) is given by 


to unity. 


a: 0° 10° 20° 30° 40° 45° 
K (a): 1.000 1.000 1.000 1.000 1.000 1.001 
M (a): 00 5.759 2.922 1.996 1.547 1.402 
a: 50° 60° 70° 80° 90° 
K (a): 1.001 1.004 1.011 1.037 1.185 
M (a): 1.290 1.141 1.089 1.084 1.185 


For a square joint of side 2a, the upper half of the adhesive layer consists of eight 
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triangular prisms in which a = 45°. With the field (16) in each prism the upper bound 


ty obtained for the average tension across the joint is given by 


ty ’ la_[mr 
—=1+ u( Reel K( 
k } 3h | (20) 
(m4 Ss 
3 OC 


This bound is 0.1 per cent above the lower bound (5) for large values of a/h. 

A polygon of n sides circumscribed about a circle of radius a consists of 2n triangles 
such as triangle OAB in Fig. 3, subtending angles a, , a2, --- , a, at the center O of the 
circle, where zs a, = m. The use of the field (16) in each triangular prism of the adhesive 


layer for a joint with such a section then provides the upper bound 


> M(a,) tan a, > K(a,) tan a, 
i oe ne (21) 


min T Sh . 


. b tan a, : tan a, 
1 1 


for the average tension across the joint. The coefficient of the term 1/3 a/h in (21), 
which is the important term for thin layers of adhesive, is less than K(a,,.), Where max 
is the numerically greatest of the angles a, , a. , --+ , a, . Thus the difference between 
the upper bound (21) and the lower bound (5) is always less than 19 per cent of the lower 
bound for large values of a/h. The term independent of a/h in (21) may become appreci- 
abie if a number of the angles a, are small and if a/h is not too large. In this case a better 
upper bound may be obtained by using the field (11) in the prisms subtending small 
angles at O. 

For elongated sections it may be advantageous to assume that regions of the adhesive 


move in a plane strain manner. As an illustration consider a joint whose section is a 
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Fic. 4. Plan for rectangular joint velocity field. 


rectangle of sides 2a, 2b, where b > a. Figure 4 shows one quarter of the rectangle. 
The velocity field is symmetrical about the lines OA, OD parallel to the sides and passing 
through the center O of the rectangle. In region O’AB the components v, , vs have the 
values (16), where the origin of the polar coordinates r,@ is taken at O’. The field in region 
O’BC is the image of the field in region O’AB with respect to the line O’B. In region 
OO’CD the plane strain motion 


v, 0, v, —Vy/h 
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is assumed. The field provides the upper bound 


ty a, laj|3 a 
y = 240.4027 + 35 E — 0.499 | (22) 
which is close to the lower bound (7). 

The upper and lower bounds obtained above determine the strength of butt joints 
whose sections are circular, rectangular, or a polygon circumscribed about a circle with 
sufficient accuracy for practical purposes. The strength of the joint is inversely pro- 
portional to the thickness of the adhesive layer for thin layers. The bounds indicate 
that, for a given cross-sectional area and a given thickness of adhesive, the circular 
joint is the strongest of all joints with a convex area of cross-section. 

Acknowledgement. ‘The writer would like to thank Dr. N. A. de Bruyne of Aero 
Research Limited, Duxford, Cambridge, for his suggestion that the analysis of Ref. 11 


be applied to butt joints. 
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Nomograms of complex hyperbolic functions. By Jorgen Rybner. Jul. Gjellerups Forlag, 
Copenhagen, 1955. 39 pp. $6.40. 
This is a very carefully prepared group of alignment charts and formulas which should be extremely 
useful in any calculations involving complex hyperbolic functions. Included are formulas for circular 
and hyperbolic functions and for four terminal networks and transmission lines. Nomograms for re- 


flection loss and phase shift and for insertion loss and phase shift are among the many alignment charts. 
Roun TRUELL 


The theory of games and linear programming. By S. Vajda. Methuen & Co., Ltd., London, 

John Wiley & Sons, Inc., New York, 1956. 106 pp. $1.75. 

The contents of the little volume are roughly indicated by the following chapter headings: An 
outline of the theory of games—Graphical representation—Algebra of the theory of games—An outline 
of linear programming—Graphical representation of linear programming (1)—Algebra of the simplex 
Degeneracy and other complications—Duality—The solution of games—Graphical represen- 


method 
The method of leading variables. The presentation is concise but 


tation of linear programming (2) 
easy to follow. In the field of linear programming, in particular, the book should be instrumental in 
dispersing the aura of mystery and complexity that this technique has acquired among many of its 
potential users because most expository material available so far either restricts itself to giving mere 


operating instructions or presupposes far too much mathematical background. 
W. PRAGER 


Frequency response. Edited by R. Oldenburger. The Macmillan Company, New York, 


1956. xii + 372 pp. $7.50. 

In December 1953 a Frequency-Response Symposium was held in New York as part of the Annual 
Meeting of the ASME. The present volume contains primarily papers presented at this Symposium. 
To round out the material, reprints of several important papers on frequency-response have been added; 
some of these were originally published in Russian and had so far not been available in English. The 
space available for this review does not permit listing of the indivudal contributions, which are grouped 
as follows, the parentheses indicating the numbers of papers in each group: Fundamentals (6)—Frequency 
Servo, Airplane and Power System Applications (4)—Process Control (4)—Transient 


Response Aids (2) 
Nonlinear Techniques (4)—Sampling Controls (2)—Statistical 


Response (3)—Optimum Controls (2) 
Methods (1). The volume is dedicated to Dr. H. Nyquist, whose portrait appears on the frontispiece. 
W. PRAGER 


Foundations of quantum theory. By Alfred Landé. Yale University Press, New Haven, 
and Geoffrey Cumberledge, Oxford University Press, London, 1955. viii + 106 pp. 
$4.00. 

‘As a whole the book is not an introduction for beginners but rather a postscript to more traditional 
books on quantum mechanics”’ (from the author’s preface). The author attempts to deduce the principles 
of quantum mechanies as a logical consequence of continuity of cause and effect; in particular, the Gibbs 
paradox of diffusion theory is removed by defining a ‘fractional degree of likeness’’ between states. 
This new approach is far removed from the traditional approach to the subject and in some places 
perhaps even is in conflict with accepted ideas. Although the book is certainly thought-provoking and 
novel, it is not likely to have much impact on the customary methods of teaching quantum mechanics 
until all possible aspects have been more thoroughly discussed. 

GORDON NEWELL 


(Continued on p. 154) 
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ON STRESS FUNCTIONS FOR ELASTOKINETICS AND THE INTEGRATION 
OF THE REPEATED WAVE EQUATION* 


BY 
KE. STERNBERG 
Brown University 
AND 
R. A. EUBANKS 


American Machine and Foundry Company 


1. Introduction. The general solution due to Galerkin [1]' of the field equations in 
classical elastostatics, was extended to elastokinetics by Iacovache [2] on the basis of a 
formal operational scheme originated by Gr. C. Moisil [3]. The method of derivation 
employed in [2], however, does not assure the completeness of the solution so obtained. 
In the present paper we supply a completeness proof for the generalized Galerkin solution 
and relate it to a generalization of Papkovich’s [4] solution of the equilibrium equations. 
With a view toward the fact that the Galerkin vector of elastokinetics, in the absence 
of body forces, satifies a repeated wave equation, we prove further that every solution 
of such an equation may be represented as the sum of two wave functions. This result 
is a generalization of a theorem due to Almansi [5] regarding the integration of the 
biharmonic equation. 

2. Completeness of the generalized Galerkin solution. ‘The displacement equations 
of motion in the linear theory of homogeneous and isotropic elastic media are given by 
eos 2. (1) 


VVa~-- 0, 
v 


Vu ; 
t hot ob 


l 
i-2 
where u is the displacement vector,” / stands for the time, F is the body force density, 
p the mass density, while »v and uw denote Poisson’s ratio and the shear modulus, re- 
spectively. Let G(P, t) be* four times continuously differentiable for P in an arbitrary 
region of space D and a < ¢ < 6. Substitution into (1) confirms that every displacement 


field of the form 


u = 5 (21 - )0%G - VV-G] (2) 
satisfies (1) provided 
noe = -—*_, (3) 
l—~y» 
*Received April 17, 1956. The results communicated in this paper were obtained in the course of 


an investigation conducted under Contract N7onr-32906 with the Office of Naval Research, Depart- 
ment of the Navy, Washington, D. C. 

‘Numbers in brackets refer to the bibliography at the end of the paper. 

“Letters in boldface designate vectors; the symbols “‘-’’ and “ X” indicate scalar and vector multi- 
plication of two vectors, respectively. \/ is the usual del-operator. 

3We shall write consistently f(P, ¢) in place of f(x, y, z, 0), P being a point with rectangular cartesian 


coordinates (x, y, z). 
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where 
- ° | a7 . 
O1=V -3>3 (¢ = 1, 2) (4) 
Cc; ol 
and 
° 2(1 — v) 2 
C; Sts 2s . C2 = e. (5) 
1 — 2 p p , 


Equations (2), (3) constitute the solution of (1) discovered by Iacovache [2], when 
the body forces are absent. If G is independent of 1, this solution becomes identical 
with Galerkin’s solution [1] of the equilibrium equations.* We now prove the following 
completeness theorem: 

Let D, with the boundary B, be a bounded (not necessarily simply connected) region of 
space. Let u(P, t) be continuous in D + B, piecewise four times continuously differentiable 
for Pin Dand —& <t < o~. If u(P, t) is a solution of (1), then there exists a function 
G(P, 1), of the same degree of smoothness, such that the representation (2), (3) holds. 


By virtue of the Stokes-Helmholtz resolution,’ there exist ¢(P, t) and A(P, f), both 
having the degree of smoothness of u(P, ¢), such that 


Quu= Vot+ V XA. (6) 


Substitution of (6) into (1) yields, with the aid of (4) and (5), 





S Dive + O3V xX A + 2F = 0. (7) 
Next define a function G(P, f) by means of 
—1 Ci 
(P = —————_ | + Vo. 8 
G(P, t) Sx(1 — ») E Vet xh], (8) 
in which ¢,(P, ¢) and A,(P, f) are the retarded potentials given by 
r (6 — R/c.) r A(Q,t — R/ ; 
(P,) = [M24 RY) a, ar, = M9, t— Bie) o, = 9) 
Jp R D R 


Here F is the distance from P to Q and the integrations are with respect to Q. Clearly, 
G(P, t) has the required smoothness. Moreover, from the theory of retarded potentials,° 





Did. = —40, OA, = —40A, (10) 
whence (8) implies 
2 Fe I ae _ 2 
OG =a» V XA Sr — nae D1 Ve (11) 


4See also Mindlin [6] and Westergaard [7]. 
‘See, for example, Phillips [8], p. 187. 
*See Phillips [8], Chapter IX. 
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By (10), (11), 


——" 1 1 
0:0:G6 = XI —» (4 DiV¢e + O2V xX A), (12) 
and (7), (12) imply (3). 
In view of (8), 
ss 
lanl 8r(1 — v)e; Vide » (13) 
while (6), (11) lead to 
2 
Qu = Vo + Al — »)OIG + 745 Diver - (14) 
2 2 
It follows from (14), (10), and (4) that 
‘ 1 (ct , . 
Quu = 2(1 — »)OiG + — (2 — VV. (15) 
2 


and (15), (13), because of (5), yield (2). This completes the proof of the theorem. 
If we define functions ¢(P, ¢) and &(P, 2) through 


o(P,t) = —V-G — r-d, W(P, t) = —40G, (16) 
in which r is the position vector of P, (2) becomes 
Quu = Vio +r-y) — 4(1 — vd. (17) 


On the other hand, we conclude from (16), with the aid of (3) and (4), that 
DetrOw=0, OF = 57—: (18) 


If ¢ and & are independent of ¢, Eqs. (17), (18) reduce to Papkovich’s [4] general solution 
of the equilibrium equations which was independently discovered by Neuber [9].’ 

3. A theorem on the integration of the repeated wave equation. According to 
(3), the generalized Galerkin vector G(P, f) satisfies an inhomogeneous repeated wave 
equation. If D is bounded while F(P, /) is continuous in D + B and piecewise twice 
continuously differentiable for P in D, —” < t <, a particular solution of (3) is 
obtainable through an iterated retarded potential. Indeed, such a solution is given by 











GIP, t) = | me, ¢ — R/ey) 4, 
D R (19) 
—] ’ FQ, t — R/e,) 
4, ( > —= — ————— inka /2 ] 
mo - GA —o I, R ar, | 


where R is again the distance from P to Q. We now turn to the homogeneous repeated 


wave equation and prove the theorem: 


7See also Mindlin [6], [10]. The body forces are assumed to vanish in [4] and [9]; they are included 
in [10). 
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Let G(P, t) be four times continuously differentiable for P in an arbitrary region D and 
a<t <b. If G(P, 0) is a solution of 


ois =0, area, (20) 
then it admits the representation 

G = G,(P, ) + GP, 0), (21) 
where 

OG; =0 (= 1,2) (22) 


and G;,(P, t) is twice continuously differentiable in D,a < t < b. 
It suffices to show that there exists a G,(P, ¢) which has the desired smoothness and 
simultaneously meets the equations 
OG, = 0, OxG — G,) = 0. (23) 
Once the existence of such a G, has been demonstrated, we define G.(P, t) by 
G(P,i) =G-—G,, (24) 


and the proof is complete. 
Equations (23), in view of (4), are equivalent to 


aG, 


OiG, = 0, a %, (25) 
where 
g(P, t) = =", 06 (26)* 
Cs — C} 
and, from (20), 
Dig = 0. (27) 
Thus, we need only exhibit a function G,(P, ¢) satisfying (25), subject to (27). 
Consider the function H(P, t) defined by 
H(P,t) =| | g(P,8) dé dn. (28) 
Clearly, G, = H conforms to the second of (25). Moreover, by (27), (28), 
——_ , 0H 
_ OH = 01 - = Og = U0, (29) 
ot ot ; 
whence 
OH = a(P) +- tB(P). (30) 


Next, consider the function G,(P, ¢) defined by 


GAP, t) = H(P, t) + A(P) + tB(P), (31) 
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where A and B are particular solutions of the Poisson equations 
V°A = —a, V°B = -8. (32) 


By virtue of (28), G, so defined still meets the second of (25), and according to (30), 
(32), it also satisfies the first of (25). Finally, G,(P, t) has the desired smoothness. The 
proof is now complete. 

The foregoing theorem is the counterpart for the repeated wave equation, of Almansi’s 
[5] theorem regarding the representability of biharmonic functions in terms of harmonic 
functions. A different kind of extension of Almansi’s theorem was given in [11]. It should 
be noted that the present theorem, in contrast to those established in [5] and [11], does 
not involve any convexity restrictions on the region D. 

As a consequence of the two theorems proved in this paper, and with reference to 
(2), (3), (4), the general solution of the equations of motion (1), in the absence of body 


forces, admits the alternative representation 


p a°G, 
at” 


oO°G;, = 0 (¢ = 1, 2). (34) 


2uu = 


—- VV-G, +G,), (33) 


This reduces the problem of elastokinetics to the determination of appropriate solutions 


of the wave equation. 
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(Continued from p. 148) 
The structural interdependence of the economy. Proceedings of an international conference 
on input-output analysis, Varenna, 27 June—-10 July, 1954. Edited by Tibor Barna. 
John Wiley & Sons, New York—A. Giuffré, Milano, 1956. ix + 429 pp. $7.50. 
“The 


twenty-one 


In addition to a prefatory note by W. Leontief, a foreword by R. Tremelloni, and a note or 


changed tasks of a faculty of economics” by G. Bruguier-Pacini, the volume contains 
ywing headings: Methods of analysis—Social accounting aspects—National 


papers grouped under the foll 
a listing of all contributions, the 


applications. Since space limitations prevent 


experience s Spe Ci i] ppl 
mly those likely to be of interest to applied mathematicians: Input-output 


if 


following partial list contains « 
Dynamic programming, by G. Morton—Industry-wide, multi- 
R i 5 2 


computations, by W. Duane Evans 

industry and economy-wide process analysis, by H. Markowitz—Input-output and the social accounts, 
by R. Stone—Aggregation problems in input-output models, by E. Malinvaud—On changes in inter- 
sectoral terms of trade, by P. N. Rasmussen—Linear expenditure systems and demand analysis: an 


application to the pattern of British demand, by R. Stone. 
W. PRAGER 


Operations research for management. Volume II: Case histories, methods, information 
handling. Edited by J. F. McCloskey and J. M. Coppinger. With an Introduction 
by the Earl of Halsbury. The Johns Hopkins Press, Baltimore, 1956. xv + 563 pp. 
$8.00. 


The reception accorded “Operations research for management” (see review in this Quarterly 14 
(1957) ...) encouraged the editors to collect a number of papers presented at the Johns Hopkins Informal 
Seminar on Operations Research during 1953-54 and 1954-55 in the present volume. While a considerable 
number of these papers were also published in such journals as the Operational Research Quarterly or 
The Journal of the Operations Research Society of America, many readers will doubtless be glad to have 
them collected in a handy volume. In addition to an introduction (“From Plato to linear programming’’) 
by the Earl of Halsbury, there are twenty-eight papers, 13 on case histories, 9 on methods, and 6 on 
information handling. Papers not previously published deal with the following subjects: Operational 
research in the British coal industry, by H. P. Privet—Operational research in underground mining, 
Design of experiments in operations 


by S. L. Cook—Utilization of training aircraft, by M. Astrachan 
-Bounding the 


research, by W. J. Youden—Operational organization, by D. B. Hertz and G. J. Feeney 
solutions of practical queueing problems by analytic methods, by G. D. Camp—Failure of complex 
equipment, by G. D. Shellard—Operational gaming in industry, by W. E. Cushen—A Monte Carlo 
model for military analysis, by R. E. Zimmerman—-How is planning possible, by C. West Churchman, 
and some of the six papers on information handling which were integrated by L. 8. Christie into a 


coherent account of work performed at M.I.T. 
W. PraGEeR 


The principles of mechanics. By Heinrich Hertz. Translated from the German by D. E. 
Jones and J. T. Walley. XLVI + 274 pp. Dover Publications, New York, 1956. 

271 pp. $3.50. 

Except for the addition of an introductory essay on Heinrich Hertz’s philosophy of science (20 


pp.) by Robert 8. Cohen, this is a reprint of the 1899 edition of this translation. 
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LINEAR PROGRAMMING AND PLASTIC LIMIT ANALYSIS OF STRUCTURES* 


BY 
W. 8. DORN 
ircraft Gas Turbine Development Department, General Electric Company 
AND 


H. J. GREENBERG 


New York University 


1. Introduction. The variational principles of Greenberg and Prager [1]' for the 
limit analysis of elastic-perfectly plastic structures require the optimization of a linear 
functional subject to linear constraints. The identification of these problems with linear 
programming problems was first made by Charnes and Greenberg [2]. In the present 
paper, the problems of limit analysis are reduced to three basic types of linear program- 
ming problems chosen so as to keep the size of the problem as small as possible. Appro- 
priate methods of solution for each of these types of problems are reviewed and are adapted 
here for the special systems encountered. Procedures for determining an initial feasible 
solution for each of these methods, including a new procedure for determining an initial 
extreme point solution for use with Lemke’s dual method [3], are discussed. The various 
linear programming solutions of the limit analysis problem are compared, and a simple 
example is carried through following each method of solution. 

The structures considered are plane pin-jointed trusses for the sake of simplicity 
of description. The mathematical formulation would be the same, however, for rigid 
frames or continuous beams and for space as well as plane structures. 

It should be emphasized that the real value of the linear programming solutions of 
limit analysis problems lies in the analysis of many variable problems arising in highly 
indeterminate structures. The methods of linear programming as outlined here lead to 
exact solutions by systematic procedures of types which are today routine for high speed 
digital computing machines. 

Two codes for the solution of linear programming problems have been written for 
the IBM 701 computer [10, 11]. Eisemann’s code [10] employs the simplex technique 
as described in Sec. 4 below. The RAND code [11] is based on the “revised” simplex 
method which has certain advantages for high speed automatic computing. This latter 
code can accommodate 200 restraint equations. A structure with approximately 200 
redundancies can be handled, therefore, using this code. 

Codes have also been reported available for the UNIVAC and the IBM 650 computer 
[12]. 

2. Limit analysis problem. Consider a plane truss with no external redundancies 
and composed of s bars and k joints. Let a number of concentrated loads bearing fixed 
ratios and lying in the plane of the structure be applied at the joints of the truss. The 
equilibrium equations may then be written’ 

*Received May 11, 1956. This research was sponsored by the United States Air Force, through 
the Air Research and Development Command. The material contained in this paper was submitted 
as part of a thesis to the Carnegie Institute of Technology in May 1955, by the first author in partial 
fulfillment of the requirements for the degree of Doctor of Philosophy. 

‘Numbers in square brackets refer to the bibliography at the end of the paper. 

*See for example, pp. 115-122 of Ref. [4]. 
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>> a;;8; = (i = 1,2,--- , 2k — 3), (2.1) 
where S, is the axial force in the jth bar, considered positive for tensile forces; p; are the 
fixed loads, and } is a positive multiplier which is allowed to increase monotonically from 
zero corresponding to proportional loading. The a,;; depend on the geometrical configura- 
tion of the truss and are essentially direction cosines of the angles between the bars and 
the coordinate axes. If s > 2k — 3 then the truss is redundant, and Eqs. (2.1) admit 
a non-trivial solution S; for \ = 0. Indeed, if there are r redundancies, then s = 2k — 


3+ f. 


If it is assumed that the bars are composed of elastic-perfectly plastic material, the 
forces are further restricted by the yield conditions 


-L,;< 8; <U Cie oy Ae ae 


where U’; (—L,;) is the yield force in tension (compression) 

By the maximum principle of limit analysis, the largest value of \ for which a solution 
S,; , to Eqs. (2.1) and (2.2) exists, is the safety factor against plastic collapse [5]. 

3. The linear programming problem of Type I. Assuming, without loss of generality, 
that the first 2k — 3 bars form a statically determinate truss, the equilibrium equations 
may be solved for the forces in the first 2k — 3 bars in terms of the forces in the re- 


dundant members. 
The maximum principle may then be phrased as follows in the standard notation 


and form of a linear programming problem. 
Maximize 


A = Pox (3.1) 


subject to 
Ps sc (7 = 1,2, +--+ , 2s), (3.2) 


where x and P; (j = 0, 1, --+ , 2s) are column vectors’ in a real, r + 1—dimensional 


vector space, V,,., . The Cartesian components of these vectors are given by 
ak 
a, ee 
’ l 
2k 5 
a 
: i 


= 2 Q ; Qj, 
’ l 
2-3 

+ D0 aiip. 


*Prime denotes transpose, i.e., a row vector, and Por, for example, denotes the inner product of Po 


and z. 


LINEAR PROGRAMMING AND PLASTIC LIMIT ANALYSIS 157 


1957 
Ps, io [6;.,°°> , 6; , 0] (j = 1,2, --- ,r) 
Pive = —P; (j = 1,2, --- ,8) (3.3) 


P? = (0, --- ,0; 1] 


xz = [8,, ppt gee Ae 


and the c,; are real scalars defined by 
\U; (7 = 1,2, ---.,8) 
c; = ¢ 
(L,. (j=84+1,-+-, 29. 


The a,; are defined by 


where 6,,, is the Kronecker delta.‘ 
Notice that. 
c, >0 (j = 1,2, --- , 2s). (3.4) 


The problem stated in (3.1) and (3.2) will be referred to as a linear programming 
problem of Type / and is naturally adapted to solution by Lemke’s Dual Method [3]. 
P,;,, = —P, , however, provides a simplification since only the first s 


The fact that = 
vectors need be carried in the calculations. 
The tableau arrangement [6] developed by Orden, Dantzig, and Hoffman is modified 
by the addition of a final row with entries —(z;., — c;,,) as shown in Table I. 
The vectors a, constitute a basis for V,,, chosen from among the P; so that z; — 
c; < 0 for all j, Pia, is the component of P, along the vector a, , and z; = 3°12} (Pla;)e, 
The last row is computed from 


—(254. — Cis.) = (; — ¢;) + (U; + L,). (3.5) 


An optimum solution has been obtained if the elements in the P, column are non- 
negative. If, however, Pia, < 0 then the negative entries in the qth row are divided 
row, and the positive entries in the qth row 
c;.,) row. The minimum of 
k, there are two cases to be 


into the corresponding entries in the z; — ¢; 
are divided into the corresponding entries in the —(z,;,, — 
these quotients is selected. If this minimum occurs for 7 = 


considered. 


Case I: Pia, < 0. Then P, replaces a, in the basis and the tableau corresponding to 
the new basis is obtained from 
r+1 , , 
Pla, P a, ” 
P,; = >>| Pla; — =i Pla; |a; + =} P. . (3.6) 
{=1 I pe Pa, 
ix@ 


Note that the a; represent the elements of the inverse of the coefficient matrix of the unknown 


forces S; , -+* , Sox_3 in the equilibrium equations. 

This permits these forces to be expressed in terms of the redundant forces in order to deal with the 
latter exclusively. However, in actual problems equilibrium equations which express each force explicitly 
in terms of the redundant forces can usually be written down directly from the structure without need 


of matrix inversion. 
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Cs « Ci Ck Cs 

a; Py P,: P, P 

Ci ay Phas Phay° id Pha: =" Pla 

Ca ag Poa, Piag° =~ Piag: les Plag 

Cra} Grit Plays Pia; 1 Pharai’ nite Pars 
| < C1 Ze — Ck By Gy 
= . a to ares 
— (2.41 — Co41) —~(Ze4e — Cork) — (Zen — Cas) 
| 

TABLE I 


Case II: Pia, > 0. Then P,,, replaces a, in the basis and the new tableau is found 


from 
: Pla Pla 
P; = Pee eo Pm i: a Paes (3.7) 
i= 2, Pia — Br per 
tF¥@ 


When the optimum solution has been obtained, the values of the forces in the members 
are determined from 
Ss; = 3 (2; oy" (Zj4. a Ci +e) + (U, rom L;) (J = 1, 2, at is , $). (3.8) 
4. Type II or simplex problem. The dual linear programming problem to the one 
defined by (3.1) and (3.2) is: 


Minimize Zo = >. pe, (4.1) 


over the variables p; subject to 


> p;P; = P, {.2) 
i=l 
and 
ep; < 0 (j = 1,2, --+ , 2s). (4.3) 


This will be referred to as a linear programming problem of Type II or Simplex 
Problem, and the dual theorem” relating it to the Type I problem states that if either 


'See, for example, Chapter VIII of Ref. [6] 
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the minimum of z) or the maximum of Pz exists and is finite, then 
Minimum z = Maximum Pz. 


The minimum problem phrased in (4.1), (4.2), (4.3) was shown to be equivalent to 
the kinematic principle of limit analysis by Charnes and Greenberg [2]. This problem 
is naturally adapted to solution by the simplex technique [6, 7]. 

Che same tableau (Table I) as used in the dual method for the Type I problem may 
be used provided the basis vectors, a; , are now chosen so that Poa; > 0 for all 7. 

{n optimum solution is at hand if z; — c; < 0 for all 7. If z, — c, > 0 then P, may 


be brought into the basis to obtain a smaller value of z, . Two cases arise: 


Case 1: 0 i; < s. The positive entries in the P, column are divided into the corre- 
sponding entries in the P, column and the minimum of these quotients is selected. If 
the minimum occurs for 7 = q, then P, replaces a, in the basis and the algorithm for 
obtaining the entries for the new tableau is given by (3.6). 

Case II: s « < 2s. The negative entries in the P, column are divided into the 


corresponding entries in the P, column, and the minimum of the absolute values of 
these quotients-is selected. If the minimum occurs for 7 = q, then P,,, replaces a, in the 
basis and the algorithm for the new tableau is (3.7). 

(gain at the optimum solution the S; are determined by (3.8). 

5. A bounded variables problem. ‘The collapse problem stated in (2.1) and (2.2) 
may be transformed into still another type of linear programming problem which 
LOLLOWS, 


Maximize al 


ipject to 





s+1 
P, = >> w/P; 5.2) 
i=1 
0 A Ww; 4 b, Gj = i. eee i oa 1), (5.3) 
here P;(j 0, 1, «++ , 8 + 1) are vectors in a space of 2k — 3 dimensions, V.4-;3, 
ind are defined as 
} a, ;L; a,;L, 
- 
>». a2;L; | a2;L,; x 
P, . . | P, = (7 = 1,32, cee , 8) 
} Aox-3,;L; | Coe-3.41y; | 
j=l J L. 
| -Pr 
| —Pa 
Fost 2 
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and w, , c; and b; are given by 


~ 1 3 
g IL + (7 2 s) 

r ) s + 1) 

A (7 L. 3 S$) 
f ) 

1 =s-+ 1) 

U; +1 1.2 ¢) 
b L ve : 

M (ij =st+1) 


where .V is an arbitrarily large positive number. 
This can be then reduced to a Type II or simplex problem by introducing non- 


negative variables y; such that 
w, +y; = b 


This results in a Type II problem of size (2k + s — 2) & (2s + 2). 

Charnes and Lemke [8] and Dantzig [9] have shown that this may be treated as a 
(2k — 3) X (s + 1) problem, i.e., the inequalities w; < b; may be suppressed. 

6. Comparison of methods. A significant difference between the simplex and dual 
methods is that in the former a certain freedom of choice may be available in the vectors 
entering the basis, while in the dual method a choice may exist in the vector leaving the 
basis. For either method, therefore, a physical criterion is desirable to guide the analyst 


in his choice. 


For the collapse problem, if at the optimum solution P; is in the basis for ) L. 2. 
+++ , 8, the jth bar yields in tension at collapse. Similarly if P; is present in the basis for 
j=st+1,s4+2,--- , 2s, then the (7 — s)th bar is yielding in compression. 


Similar interpretations are available in the bounded variables formulation (Sec. 5 

These facts allow the experienced analyst to use his judgment and intuition in making 
choices regarding basis vectors 

A more quantitative comparison of the methods of solution arises from a consideration 
of the number of arithmetical operations necessary to achieve a solution. 

As a measure of the arithmetical operations, the number of multiplications to be 
performed per iteration will be used. For Type I and II formulations this measure is 
(r + 2) (s + 1), while for a bounded variables problem the number of multiplications 
is (s — r + 1) (s + 2). The preference for formulation on the basis of the number of 
arithmetical operations depends, therefore, on the relationship between s, the number 
of bars, and r, the number of redundancies. 

7. Initial solutions to the programming problems. In this section methods for 
obtaining initial feasible solutions to Type I and Type II linear programming problems 
are presented. For bounded variables problems it is sufficient to find a basis for the 
vector space of the equalities (5.2). The method developed for Type II problems produces 
such a basis. 

A. Type I problems. If any point x satisfying (3.2) can be found then a simple change 


of variables will translate this point to the origin. Starting from the origin, the following 
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method then produces an extreme point solution (or initial feasible solution) by use of 
the dual method applied to a modified problem. 
lor the collapse problem, c; > 0 and this implies that « = 0 is a solution to (3.2). 
Consider now the modified" problem to maximize Pjx subject to 


Pie Ss ¢ () 1,2, --- , 2s) (7.1) 
ae'x <0 (§ = 1,2, +--+ .¢ +3) (7.2) 
c, > 0, (7.3) 
where e, is a unit veetor in V,,, with +1 as the 7th component and ga; is either equal to 
| on | and is chosen according to the criteria described below. 
Notice that the origin is an extreme point of the modified problem regardless of the 
choice of sign for ¢, . The basis vectors associated with this point are e; for? = 1,2, --- , 
| 
Phe set of points x satisfying (7.1) is designated by A. The set satisfying both (7.1) 
and (7.2) is A*.*Note that (7.2) is just a restriction to some orthant of V,,, once the 
signs of o, have been chosen. In order that the solution to the modified problem coincide 
vith the original problem, it is necessary and sufficient that the ¢, be chosen in (7.2) 
so that the point n A for which Pir takes on its maximum also is contained in A*, 
1. the correct orthant of V, , must be chosen. 
The advantage of the modified problem is, of course, that « = 0 is an extreme point 
ition. Starting from this solution and using the dual method, a value of o; for some 
2 is chosen at each iteration, and the corresponding e,; leaves the basis in favor of some P, . 
he procedure for accomplishing this is as follows. Let the basis at some stage be 
k's e,, Py... +++, P.., where the signs of ¢,,,; , +++ ,o,., have already been properly 
chose 
Define r 4 vectors a; (j = 1, 2, +--+ ,r + 1) such that 
i @ j - 
Ca } (7 1) 
lO if i#J 
lol z un 
? r & Fs ] — 
Pla (7.5) 
0 if t~# J 
fon t+ 1,- rt+il. 
Chen 
aad 
P, = > (Plaive; + >> (PladP; 
! i t l 
and if .v, is the extreme point of A* corresponding to this basis then by definition 


e'x, 0 (7 a ee - 
; (7.6) 


Che authors are indebted to Dr. C. E. Lemke for suggesting this modification. 
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< i — oe, ee ae 


Pot = > (Pha;)(e:xr%o) — 8 5 (Pha,)(e'a,) + (Pha,)(e2#) 
—— s ’ , Y d 


t= ; 
r+1 
+- 2. (Pla;)(P!x) — 6 p (Pja;)(Pia,)- 
i t l t=t+1 
Now from (7.4) and (7.5 
ela 0 (2 l q-—-l,q+! 8 


and using (7.6 


The last sum is P42, so 


Piz PsP +- (Pha,)(e,z). 


0; (ii) Pha, > 0; (ii) Pha, = (), 


Three cases arise: (i) Pa 
. It follows that any # yielding a larger value 


For case (i) if ef ¢ > O then Pjt < Pox 
of the functional than x, cannot lie in the half-space e(4¢ > 0, but must satisfy 


+e’t < 0. 


Thus Ca is removed from the basis and we pick o; = +1. 
For case (ii) if e# < 0 then again Pj# < Pox, . Thus similarly it is necessary that 
—e,t < 0. 


is removed from the basis but Un is chosen to be —1. 
Pix, and the choice of o, is deferred for the present. 
remaining in the basis 


Again ¢ 
Finally for case (iii), P62 
If at some stage of the computations, case (iii) holds for all e; 

is arbitrary. 

7) are used to proceed to a new tableau. 

are eliminated from the basis and the 


then the choice of oc; for those « 
The usual algorithms (3.6) and (3. 
In this way in r + 1 iterations all of the e, 


appropriate values of c; equal to +1 or —1 are chosen. 
B. Type II problems. Consider the following modification’ of the Type II problem 


2), (4.3): To minimize 


formulated in (4.1), (4.2 
> pc; + DviM 


~] 
~J] 


subject to 


2 r 1 
Po = > oP, + Dd ve. (7.8) 


7This modification was first suggested by Dantzig [7], see footnote to p. 340. This is the technique 
employed in the 701 codes [10, 11] to obtain initial solutions, 
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o>0 (j=1,2,--: , 2s) 
7,20 (2 L. 2 «+ .f + 3), 


(7.9) 


where e; is a unit vector with plus or minus one as the ith componeit depending on 
whether the 7th component of P, is positive or negative, and where . is an arbitrarily 
large positive number. 

The solution to this problem is identical with the one phrased in (4.1), (4.2), (4.3) 
since the minimum will occur for y; = 0 for all 2. 

A basis for this problem, however, is readily available. Indeed, the basis a, may be 
taken to be 


a; = e, (7 ee Ses 2 oo i 
It follows from the definitions of e; and a; that 
a e; (2 Lie ge 2 


The entries of the tableau (Table I) are, therefore, easily computed. 

The vectors e, need not be carried in the tableau since if an e, leaves the basis it 
cannot return because it carries a large positive weight, 4/7. It will require exactly r + | 
iterations in order to obtain a basis comprised entirely of vectors chosen from among 
the P; 

8. Examples. To illustrate and compare the methods of solution of linear pro- 
gramming problems described in the previous sections as applied to structural collapse 
problems, a simple example and its solution by these methods are presented here. 

Consider the once redundant truss in lig. | loaded with a single concentrated force 





N 
\ 











G 3 
© | 
. ™ 
a Sa 
Rn 
Fig. 1 
as shown. The members are numbered as indicated, and the equilibrium equations may 
be written 
—S, — £S,; = Ab, 
S. + $S, = 0, 
S. + 4S, = 0, 
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The fully plastic forces in tension and compression are taken to be the same and equal 


to N, . The yield criteria (2.1) become then 


S| = (2 Oe ae |e 
Now letting x, = S,/V, and x, = Ab/N, and solving the equilibrium equations for 
S,, S,, °°: , 8, in terms of these, the yield criteria may be written 
— $501 + 4% | <1 
—§7| <1 
< | 
+7 T <. l 
5 4 
l < ‘, 
< | 


The safety factor against collapse is the largest value of x..V,/b consistent with the above 


inequalities 


Define 13 vectors a two-dimensional! space as 


P! P! Ae 3 ; P! P! -, 0), 


9° 4 


P! P! —3,0), Pi=-P! ‘ 1), 


42 
| Hi 
va) 
D 
1(0,4/5 
a 1 es eee en Ss 
A(-1,0) Xy 
(O1- 4/5) 
iy, 


Fic. 2. 
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The problem in vector notation is then to maximize xr, = Px subject to 
Piz s 1 (j = 1,2, --- , 12). (8.1) 
The equation Pf. | defines a line in two-space for each 7. The vector P; is normal 


to this line and points into the half-space for which the corresponding inequality is 
violated. All of the lines defined by equality in (8.1) are shown in Fig. 2. The set of points 
A for which all of (8.1) are satisfied in the parallelogram ABCD, and the maximum value 


of x, is obtained at the point D (1, 8/5). The safety factor is then 
8 N, 
{= == 
ov b 


\n analytic solution is obtained by the dual method. Let 
l 0 
li 
0 l 


be the initial basis for the space. This is not an extreme point solution of A but is used here 
to find such an extreme point solution. The initial tableau appears in Table ITA. 

Since Pha, > 0, ¢ is chosen to be —1. This multiplies the e, row by —1. Now choose 
the minimum over / of 


2(P;) ; 
p’ for Pla, < 0, 
QA» 
—2AP 26) , 
P’ for Pla, > 0. 
Qs» 
This occurs for ; 5 and is 4/5. Since Pfa, > 0, P,, replaces e, in the basis. 
Phe result is Table IIB. The process is iterated and o, is chosen to be —1 after 
which P, replaces « Table IIC represents the final solution. 


The maximum value of ., , 8/5, appears in the P, column and the z(?;) row. Thus 
8, /5b. The presence of P,, and P,, in the basis indicates that bar 6 yields in tension 


and bar 5 in compression, i.e., S N,,, Sy = N,. Because of the normalization used 
here, Eq. (3.8) for the SS, is 
S Ml(z; — ¢;) — (245 — C)4)]N, 
Chu 
S, 3N,/3; S, —4N,/5, S —3N,/5, S, = 4N,/5. 


lhe dual to this problem is a simplex problem and the same initial tableau (Table 
[1A) may be used if the vectors e, , e, are given large positive weights, 17. Now many 
ectors, i.e. those for which 2(P?,;) > 0, may enter the basis. P, is chosen on an intuitive 
basis and thus ¢, leaves the basis. The resulting tableau is shown in Table IITA. P, then 
replaces e, and ‘Table IIIB results. Finally P,, replaces P, and the final tableau will be 
identical with Table IIC. 

Notice that the solution in Table IILA corresponds to the point (1, 44) and the 
solution in Table IITB to the point G(1, 9/5) in Fig. 2. 

Finally, the problem can be formulated as a bounded variables problem of size 
BD We ai. 
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( ] l ] ] l l 
| 
} 
| J a P P P P; P; P P 
| 
0 0 3/6 | —4/5 | -3/5 | —4/5 
\ 
Z(P 0 l ] —|] l —| 
—Z (I +] l l +1 (+) +1 
C; | l l | l l 
| 
| 
| 
| 1 P P, P P; P P 
— 0 1/5 0 1/5 3/5 0 0 l 
B 
l P 1/5 3/5 U 0 1/5 | 0 
| a(P 1/5 2/5 l —] 1/5 2 C) 
| 
—Z(P 8/5 +] 9/5 0 +] 
C: I I I I 
| a P P P P; P, P; P, 
1 Pr 1/5 0 1/5 3/5 0 0 1 
( 
l P 1/5 3/5 0 0 1/5 I 0 
Z(P 8/O 2/5 9/d 8/5 1/5 2 0 
—Z P 5 4¢ 8/5 1/5 +2/5 9/5 0 +2 


TABLE II. 
Dual method solution 
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Surveys in mechanies. The G. I. Taylor 70th Anniversary Volume. Kdited by ¢ 
175 


Batchelor and R. M. Davies. Cambridge University Press, New York, 1956. vii + 


pp. $9.50. 

This book has been written to commemorate the 70th birthday of Sir Geoffrey Taylor. It contains 
a biographical note by R. V. Southwell and ten articles in the form of surveys of the present position 
in some of the fields of mechanics in which Sir Geoffrey has been active and to which he has made im- 
The Mechanics of Quasi-Static Plastic Deformation in Metals” by R. Hill 
he phase of the development of the subject from 1950 to the time of writing 
2. “Dislocations in 


portant contributions. 1. 


the author concentrates on 1 
The Principle of Maximum Plastic Work and its Applications”. 


and emphasizes 
a general survey of the achievements of the study of dislocations. 


Crystalline Solids” by N. I Mott 
3. “Stress Waves in Solids” by R. M. Davies; a review 
dispersion of elastic waves in a circular cylinder and of visco-elastic a 
Fuilds” by H. B. Squire; a review of the major theoretical developments of the motion of an ideal fluid 
which is assumed to rotate as solid body in its undisturbed state. 5 The Mechanics of Drops and 
Bubbles” by W. R. Lane and H. L. Green; the production of liquid drops and their motion in a gaseous 
s bubbles and their motion in a liquid medium. 6. ‘‘Wave Generation 


ol experiment il techniques and accounts of 
nd plastic waves. 4. “Rotating 


vn 


medium and the production of 
by Wind” by F. Ursell; a critical discussion of certain aspects of the subject: instability calculations 


semi-empirical theories of “Viscosity 
Effects in Sound Waves of Finite Amplitude” by M. J. Lighthill; an account of the influence of thermo- 
viscosity, heat conduction, relaxation) and of nonlinearity effects 


mphasis on the conflict of these two influences. 8. ‘Turbulent Diffusion”’ 


wave generation and wave forecasting and sea roughness. 7. 


dynamically inversible processes 
on sound propagation with « 
by G. K. Batchelor and A. A Townsend; a review of the present state 
essential ideas rather than the empirical working rules. 9. ‘‘Atmospheri 
f turbulent air flow over the 


of basic knowledge of turbulent 
diffusion with emphasis on the 
Turbulence” by T. H. Ellisor 
ground. 10. “The Mechanics of Sailing Ships and Yachts 


pal features ol yachts and high speed sailing craft 
G. MorGan 


1 discussion of the empirical knowledge o 


’ by K. S. M. Davidson; an account of the 


development of the princi 


Journal of fluid mechanics. Editor: G. K. Batchelor. Assistant Iditors: T. B. Benjamin, 
I. Proudman. Associate Editors: G. F. Carrier, W. C. Griffith, M. J. Lighthill. 
Taylor & Francis, Ltd London, and Academic Press Inc., New York. Subscription 
price $16.50 per volume of approximately 600 pp. 

levoted to theoretical and experimental investigations of all aspects of 


The new journal is to be ¢ 
volume of approximately 600 pages. The first 


mechanics of fluids, and is to be issued in six parts pe! 
part (May 1956) contains papers by M. D. Van Dyke, P. 
L. c. Woods, V. Blackm in, J Crease, M. B. Glanert, ind N tott 


G. Saffman and J. 8. Turner, M. J. Lighthill, 


Integraltafeln zur Quantenchemie. By H. Preuss. Springer-Verlag, Berlin, Gottingen, 
Heidelberg, 1956. 162 pp. $9.35. 


This book is concerned with the calculation of wave functions and allied quantities arising in the 


binding and molecular physics. The first part, of three parts, is con- 


quantum mechanics of chemi 
The second part of the book is 


cerned with the formulation of problems and the methods of solution 
concerned with methods of handling the integrals involved and the third part contains over 100 pages 
of tables of functions used in evaluating the integrals of parts I and IT. This is a handbook for the person 
doing serious calculation in quantum chemistry. 
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ASYMPTOTIC SOLUTIONS OF A CLASS OF ELASTIC SHELLS 
OF REVOLUTION WITH VARIABLE THICKNESS* 


BY 
C. NEVIN DeSILVA anp P. M. NAGHDI 


University of Michigan 


1. Introduction. Since the pioneer work of H. Reissner [1] and Meissner [2, 3] on 
the small axisymmetric deformation of thin elastic shells of revolution, numerous in- 
vestigations on this subject have dealt with shells of uniform thickness, but comparatively 
little attention has been given to shells of non-uniform thickness. Meissner in [3], follow- 
ing an analysis of shells of revolution of variable thickness (to which further reference 
will be made presently), treated the case of conical shells of linearly varying thickness 
in detail. Other contributions on shells of revolution of variable thickness have been 
made by Ekstrém [4] and Spotts [5] for spherical shells, and by E. Reissner [(6] for shallow 
shells of revolution. In a recent paper by the present authors [7], the differential equations 
of shells of revolution with small axisymmetric displacements, as given by E. Reissner 
[6], were combined into a single complex differential equation, a solution of which, valid 
at a turning point’ of the differential equation, may be obtained by a more recent method 
of asymptotic integration due to Langer [8]. These results were subsequently applied to 
ellipsoidal shells of uniform thickness [9], yielding a solution valid at the apex of the shell 
vhere a regular singularity occurs in the differential equation. 

In the present paper, with reference to the differential equations given in [6 
a class of shells of revolution of variable thickness is further examined. First, the impli- 
cations of the thickness variation specified in [7] which led to the complex differential 


equation in a form amenable to treatment by Langer’s method of asymptotic integration 





and [7], 


ure studied in detail and then asymptotic solutions are given which are valid at turning 
points ol the differential equation. 
2. The basic equations of shells of revolution. With the use of cylindrical co- 


ordinates r, 6, 2, the parametric equation of the middle surface of the shell may be 
vritten as 
r r(é), z = 2(g). (2.1) 


Denoting by ¢ the inclination of the tangent to the meridian of the shell, then 


yr’ a Cos ¢, 2’ = asin do, (2.2) 
, a/¢’, r, = r/sin d, ; (2.3) 
v he re 
2 9\231/2 ¢ 
a io’) +(e’) Ts (2.4) 


*Received June 19, 1956. The results presented in this paper were obtained in the course of research 
sponsored by the Army Office of Ordnance Research under Contract DA-20-018-OR D-12099 with the 


University of Michigan 
Such points are defined here as ones at which ¥? (see equation (2.6)) vanishes to some degree 


2 and/or A has a simple or a double pole. It should be mentioned that this definition of turning 
illed transition point) includes the usual definition given by Langer where y is re- 


I 


point (sometimes « 
stricted to be a positive integer. 
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r, and r, are the principal radii of curvature of the middle surface, and prime denotes 
differentiation with respect to &. 

The appropriate expressions for the stress resultants NV, 
M, and M nd the transverse shear stress resultant Q are given in [6] as well as [7] 


A 
and will not 


and N, , the stress couples 


be repeated here. We recall, however, that as in [6], it is convenient to 


express V; and @ in terms of “horizontal” and “vertical” stress resultants, H and V, 


given by 


al r’H + i al) = - 2/H — r’V. 2.0 


ween the stress strain relations and the differential equations 


By proper elimination bet 
of equilibrium and compatibility, E. Reissner in [6] deduced two coupled second-order 
differential equations governing the small axisymmetric deformation of shells of revolu- 
tion which, while they differ only shghtly from the previous formulation of the theory 
due to H. Reissner and Meissner, contain certain advantages. It was shown in [7], that 
the differential equations of shells of revolution, as given in [6], may be combined into 
the following mplex d fferential equation: 
\ 1/2 
We’ + [eww + AlE)|W *. =) f(é)[F + 7kG], 2.6) 
age JS) | 

provided /:, given by 


is a consta 


The various quantities appearing in (2.6) and (2.7) are defined by 
: Swe la oe mrH 
Ww =| ") @+iky, yv="> 
E 7 = 
vr )(%) 
s(t k+ >} — t) 
: ( ' “aN s 
a aan (7 a)!’ I | ay | _ (“) _ 3 Bi /o)’ h’ 7 3 (3) a 
i 2 (r/a t | (r/a r | 2 (r/a) h 4 \h 2 Ih : 
8 
"f(€) = < “ ; m- IA@1l—-—y 


h f(r? Jax)! rin 
d os a ‘we oh sok ae 7 
A \ (7/a) r hj 
be (4 ao la 
Se on ‘7 _ j 9) Bea 
§=2 f + 2 le ), 
E Ag h r h (r/a) h 
where 6 is the negative change in ¢ due to deformation; h, is the value of the thickness 
h at some reference section; F and G, as in [7], are functions of the load intensity; and 
E and v are Young’s modulus and Poisson’s ratio, respectively. 
As pointed out in [7], the condition imposed for the validity of (2.6), namely that k 


or equivalently (vA — 6/2) is a constant (say — K/hp) results in the following differential 


equation for the thickness h: 
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(nr 4 (2)'n — | (Ew + (Ey = (*\ks@) (2.9) 


the solution of which, since 


j r m 

| sa) * ag = 2, 
« a u ho 
may be written as 


( . - 
} 


) 
h : ry Ky | r "Po dt +e, + 2 | ioe déf, (2.10) 


where c, and c, are constants of integration, and K, = Km/u*h, . It should be mentioned 
here that the form of thickness variation specified by (2.10) was first obtained by Meissner 
[3]. 
In order to isolate the terms in the coefficient functions of W which involve derivatives 
of h, we modify (2.6) into a new normal form as follows. From (2.9), 
: E + ] » «Set. 2, toe. 2,2 (2.11) 
2. a (r/a) h 2 2 (r/a) 2 rh 


which expression occurs in A. Also, the function ¥ in (2.8) may be written as 





2 _ 1 ho -¥ Ir'/a)’ path 919) 
y = k Wp 4% {Cle +3 a \ (2.12) 

Substitution of (2.11) and (2.12) into (2.6) after considerable manipulation results in 
W" + [Pui + (Ao + AD)IW = R, (2.13) 

where 
2_( 3 Ka\(he) ” 

WV. = L — —_ f< : 
y’ (i 58 ae, 2.14) 
_ _3(rV_1(a’\’ 1 («“) (1 + ») (r’/a)’ . 
oan (“) ~ (2) tT o\a 2 (r/a) ’ (2.15) 
3h’ (" l m') — 
A, 2” h\r 2v h (2.16) 


and FR denotes the right hand side of (2.6). 

We choose £ = ¢ (i.e., 7; = a), and proceed to examine the character of h as given by 
2.10) in the neighborhood of @¢ = 0, and the corresponding behavior of the coefficient 
functions of W in (2.13), with the aim of obtaining solutions of (2.13) valid at turning 
points of the differential equation. To this end let the principal radii of curvature be 
of the form 


= ¢" Dp”; 
— (2.17) 
r= ¢" dat’, 


n=O 


r 


~ 
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where the power series representing non-vanishing analytic functions are uniformly 
convergent in the real interval J,:0 < @ < ¢*. Substituting (2.17) into the expression 


for r’ involving both r, and r, (by (2.2) and (2.3)), and equating the coefficients of like 


powers in ¢ lead to two cases. In one case b, = —1 (and b, 0, 1) and in the other 
b, = b, = b # —1. These cases are considered in detail in Secs. 3 and 4. 

3. Case A: Behavior of shell thickness specified by (2.10) and solution of Eq. 
(2.13), when b, —1. In this case two possibilities occur: (1) }, —b, = 1 and 
4g, = 0, and (2) b, = —1, b, = 0, and q, = po. 


(1) For the first possibility, it follows that 


r= To’ (3.1) 


n=0 


and by (2.10 
h= > he’, (3.2 


where 


do 
qi =, ga=H, q2 = Qa — z : ete (3.3a) 
‘ l 0 
h=hO =a, h=0, h-= (on, oe ), ete.  (3.3b) 
do 2 ~ do 
We note that lim,.., h’ = lim, (2h. ¢) = 0 and that A in (3.2) is a non-vanishing analytic 


function. By (3.1) and (3.2), 


‘ > (n + 2)q*..9" 
r a=0 
inl co ; 
> ae" 
esi (3.4) 
Lt Dhaesd 
h . Q ’ 


and from (2.8) 


yo = m (P:)(Pe), f : b Ds “| (3.5) 


i. 
xan? 


Substitution of (3.4) and (3.5) into (2.14), (2.15) and (2.16) results in 


— 2 «as 

Vv; og VF’, 7") = § — =-141— 
2 m seit 
(3.6) 


= lp, .-1 
A ——g? — se" 4+ at, 
2p 
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where V¥ and A*¥ are analytic in J, . Similarly since 


lim A, = lim 3y (2 “s hays (3.7) 
6-0 6-0 hy \qo v ho 
A, is also analytic in /@. 

Following Langer [8], in view of the character of the coefficient functions Wf , A» , 
and A, (as given by (3.6) and (3.7)), a solution of the homogeneous differential Eq. 
(2.13) asymptotic with respect to n” (k —37/2 K,/m) to the true solution is* given by 


W = oo wt? @' [AJ o/(9) + BI-2/s(n)] (3.8a) 
where 
n = iu, > = | t*’*w*(t) dt (3.8b) 


and J and Y (to be introduced later) denote Bessel functions of the first and second 


kind respectively 

Evidently neither h, nor its derivative, influences the order of the Bessel functions in 
(3.8a). Hence this solution is also valid when h is uniform, provided k is approximated 
by a constant. It should also be mentioned that the behavior of h as given by (3.2) and 
(3.3) is similar to that discussed by Meissner [3], where r was assumed to have the form 

q cos} 

2) With b 0 and b, = 1, again r is given by (3.1) and (3.3a), and the function 


h is of the form (3.2) but the coefficients h, become 


ho=h0) =ag, h=»%+e,™, ete. (3.9) 


Jo Jo 
It is clear that h, given by (3.2) and (3.9), as well as A’, are non-vanishing analytic 


functions in [¢ 
While uy’ is again given by the first of (3.5), the function f, by (2.8), becomes 


Pia 
-n=0 Po JS (3.10) 


> ee 


n=0 Yo 


f=¢ 


and since g, ~ 0, the expressions which correspond to (3.4) and (3.6), are 


. > (n + 1)q¥*.,¢" 


Ce a 0 A 
r x r , 
> a 
oe (3.11) 
y! » > (n + I)h,.10" 
h — 


24 statement of Langer’s theorem may be found in Sec. 5 of Ref. [10]. 
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| sx, 
Vi= evr, wo) (i — 54%), 
2 m 
(3.12) 
o hy q l hi) 
10 = 8, (0 
aa 2° h Jo 2v ho/ 


and the function A, (not recorded) may be shown to be analytic in J¢. Again, an 





asymptotic solution of the homogeneous equation (2.13) with respect to yu [(k — (3/2) 2 
(K,,/m)] is given by 
W = oo yt'?@'27 AJ, a(n) + BU, (nm) ], (3.13a) 
where 
> ied 
n=07ub, b=] Uwe) dt. (3.13b) 


70 


It is apparent that the remark which follows (3.8) concerning shells of uniform 
thickness is equally applicable to the above solution. Before closing this section, we 
illustrate the foregoing results [Eqs. (3.9) to (3.13)] for a toroidal shell, where the principal 
radii of curvature of the middle surface may be specified by 


r; = R 9 


(3.14) 

1 P R 

m™=¢o al |——]+—-¢o¢], 

sin ¢, a 

where R is the radius of the circular cross section and a is the distance from the center 
of the cross section to the axis of revolution. It follows from (2.17), that 

) R, ,, = 0 for a 4; ~— 
I I _ (3.15) 


Jo = a, qa =R, etc. 


and by (3.5) and (3.10), 


(£2) 

Lu ) ‘. m, — 
P ‘ . 16) 

f = $4 - 5 E 43 sin |} 

(SIN @ a J 


Hence the function ¥* and ® in solution (3.13a) become 


a ? (Hef ¢ | R. | 
VF =(|k- .—— | — Kh -| ‘ Po 
i () 2 mt h \sin d i er J (3.17) 
/ 3 K,' . we 4 ( h\( f ) R. | 1/2 
>= [Kk - ben + —sin ? dt, 
¢ (7 9 ot | t h, \sin 7 1-4 7 sin f J t 


#0 “ 


where / is given by (3.2) and (3.3). 
For toroidal shells of uniform thickness the above solution reduces to that given by 


Clark [11], in which case h = ho = K,R/v conforms to (2.10) with c, = c, = 0, and 


‘ | » 2°]1/2 
/ tizt+]1-— (~) ‘ 
A B 
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4. Case B: Behavior of shell thickness specified by (2.10) and solution of Eq. (2.12) 
when b, = b, + 1. In this case, with b, = b, = b # —1, it follows from (2.2), (2.3), 


and (2.17) that 


Piit+db, B=24+d, _ ete. (4.1) 
Jo qi 
and by (2.10 
Cc . : -Abions« > b fe l 7) be 7 4 
} —< Oe ee lead ee Ee 0 SF 30) (4, 
F’,(¢) + e,qQ6¢ F (od) + Kyp b+2b4+2—-rb4) 1 t 2) 
d 
( in , b + 1 b+1 
/ © ( ) “ld - - ess | wee eer ( ) 
l Jol b)vp F.¢) + Kip tare —yb a? F(@) 1.3 
lhe functions / |, 2, 3, 4, 5) in (4.2) and (4.3) are defined by 
[ @ | a do - =. fin? ; 
as ; (4.4) 
; ; 7 v(1l + b) Pi a] q 
(0 ' — - * -, te. 
I l, fos 9 | tne eS et 
F.(¢) = [q%@ tT’, F.0) = 1, (4.5) 
F mre o ? [ oz d F,(0) = 1 (4.6 
(od) . pe — P (2 "(Q) = 
o Po eee es ee ee ? r az dd, 3\U) ’ \+.0) 
F.(¢) = Fle) + —*—— Fi) 
1 + b)p - 
‘ (4.7) 
@G; +. b)p > 
l+p be oe Ls + ¢0(0(1), 
(] + b)p 
f Q I @ T h e 9 I 3\D) 
- (4.8) 


34 5p {(2 + bd) — e(1 + b) E 5 + 2b ; a] q\ 
— <4 —_— ————— sss a —_ _ um ) = 
2+6(\(3+ 6 —vnl+bLp3+ 56 ali ti qo +? q % 


+ ¢'0(1), 


where 0(1) denotes a bounded function. 

With K, + 0, examination of (4.2) reveals that h is bounded at ¢ = 0 if b > —2 and 
either c, = 0 or (1 + b) v > O. Likewise with K, + 0, h’ (0) is bounded provided b > —1 
and either c, = 0 or (1 + b)v > 1. Thus, the requirement that both h and h’ be bounded 
at d 0 when K, ¥ 0 gives rise to the following restrictions’. 


b = > .9¢ 
> —1, v>O0O, K, ~0 (4.9a) 


c, = 0 or (1+ by > 1. (4.9b) 
’Although in the linear theory of elasticity, Poisson’s ratio » may have the range —1 < » < 1/2, 
as is usual on physical grounds, the negative values of » are ruled out here. 
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It may be noted that, in addition to conditions (4.9), if in (4.4) we set f,, = 0, then 
h’ (0) = 0. 
If only condition (4.9a) is imposed on (4.2) and c, ¥ 0, then h is continuous in i, 


and h’ at @ = 0 reduces to 


lim h’ lim = fir + vergo(l + be?" ™"?. (4.10) 

It is clear by (4.9b) that, if (1 + b)vy <1, then A’(0) will not be bounded and, consequently, 

a drastic change of h will take place in the neighborhood of ¢ = 0. In the case of ellipsoidal 

shells, for example, this situation leads to the occurrence of a cusp-like variation in h 
about @ = 0. 

With the restriction imposed by (4.9a), i.e, b > —1, » > O, it follows from (2.8). 


(2.17), and (4.1), that 


Peg 


tag 


We now proceed to examine the character of A, and A, given by (2.15) and (2.16) in 
the neighborhood of ¢ = 0. For this purpose, we record the following expression: 





3. h’r’ 3 (! a. | by” e195 
o"her )=(6«4NA OD” h, ® 
L+ ¢70(1) — 9g" '"'l + 01 
j 
3r(1 + b ' . 
=e 9 P d0( 1) 
K,1+ 5b Pp 
‘ L 0 (4.12 
P . AL bh? | o0( 1 1.] 
C19 K, p 1+ b 
— pO) os 0) 
h ? [ ? | h j = + § — p(l + D 
| i+ be 
i o0(1 Bess T eo d {(? q 
J j Hi P q 
. gq, 1 + (1 + 0b): 9,, ¢ Jog: 1 + (1 + bpp ae 
’ 2v > G,( 
“—— 3 7 1 + bp } v h q | + Q J + ?—) 
and write A, as 
A; B 
\ +—!4 at, (4.13 
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where 
A, = -i (40° + 8b + 3), 
(4.14) 
_ 3 (2: a) b p, 
B, = 9 (1 + b) . . 4 2 pa 


and G, and A* are analytic in J¢. 
Denoting by A the sum of A, and A, , then with the aid of (2.16), (4.14), and (4.13), 
we have 
At , BY 


A=75 - 4.15 
re aa r 4. (4.15) 


where A, is continuous in /¢, and 


‘x , 3 2 C190 l+bh)y Odo l+h)e 





ry - kK, | “4 b Po 1 b 
Bt = B, + S01 + WSs + 2 _—- ¢ 
2 , hy ll 2 ran . 
/ i +b (1 + b)p (4.16) 
7 C1Go (140) A, mM i+ os 
“tg [ Tho fu 2+ b— (1 + dp® 


l +b CiJo cise [(@ qi gait a + be) 
re fir hho ? Py, . qo (1+ bpp 


9, 190 i LO + bY sy, \ 
- hy (1 + b)p y 


If now condition (4.9b) is also imposed, two possibilities will result which will be 
considered separately: (1) c, = O and (2) (1 + b)y > 1. 


1) With ¢ 0, the functions A* and B* in (4.16) simplify as follows: 
2% 1 
Ms ; 1+ OK, Do ba - 
3% a 4 ql b . . .17) 
ij I 9' 1+ rf + [ h2+b6-(1+ by? (4.17 


Since B* , as given by (4.17), makes no contribution’ to the asymptotic solution of 
(2.13), it follows that, in the neighborhood of @ = 0, the behavior of A is dominated by 
1,,¢. Also the function ¥; by (2.14) and (4.11) may be written as 
a 3 .K, ; 
Vv; o VT (>), vr(0) =k—- 51 (4.18) 
2 m 
Hence, with (4.15), (4.17), and (4.18), an asymptltic solution of (2.13) with respect to 
u-{] 37K,,(2m)] is given by 


ae ee ee Assen add + BY eacn-dd), (4.19a) 
This may readily be seen if one introduces the transformation ¢ = s?/4, W = s'/2 U into (2.13), 


in which case the transformed A reads as (—3/4 + 4A*)/s? + B¥ + s?/4 A, where the last two terms 
are continuous in / 
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where 


n=i ub, b= |] CUA dt. (4.19b) 
(2) Here Eqs. (4.11) and (4.16), which have been deduced under condition (4.9a 
i.e., b > —1 and » > 0, are to be further restricted so that (1 + b)v > 1. 


Since for (1 + b)v > +2, A, is continuous in /¢, solution (4.19) is again valid in the 
= A, and hence, (4.19 


neighborhood of ¢ = 0. Also, when 0 < (1 + b)v < 2, limg.o A¥ = 
is a valid asymptotic form of W in the neighborhood of ¢ = 0. 
As an illustration of the above results consider an ellipsoidal shell, the middle surface 


of which is specified by r°/a” + 2°/e° 1. The principal radii of curvature may then 


be written as 


r, = a(c/a) fl — ¢ sin’ ¢] ae (4.20a) 
ap a(c/a) [1 — e¢ sin’ ¢] sas 
where 
e = 1 — (a/c)’. (4.20b) 


Comparison of (4.20) and (2.17) vields 


b=0 
Do = ale a : n= Q, Po > (3a/2)(c/a ‘@? etc., 1.21) 
GJ = alc/a gq, = O, q a/2)(c/a) ‘« et ¢ 
and by (4.11) and (2.14 
m(a/ho)(e/a)~ 
u TTi\ ¢ / } ‘ 
(4.22a 
1 — € sin’ @) 
9 , 
Stee wen ) ; ; . 
vr = VT (7 = 1 (- } (l —e sin @d) . (4.22b 
2 m f) 
Thus by 1.19)"we finally have 
WW A pP 1J,(m) + BY,(n)] (4.23) 
where 7 and 9 are evaluated by 1. 19b and }.22bD). 

It may be noted that for ellipsoidal shells of uniform thickness, solution (4.23 
reduces to that given previously in [9], where k [defined by (2.7)] is approximated to a 
constant. In the case of spherical shells of uniform thickness, (¢/a = 1), however, no 
such approximation is necessary, since / i (v/u) + [1 vip) |" and h h 
(K,/v)a. 


\ 


5. Reduction to the theory of shallow shells. We conclude the present paper with 
the special forms of the several solutions given in Secs. 3 and 4 (cases A and B) approp- 
priate in the theory of shallow shells of revolution, where attention is confined to small 
values of the meridional coordinate ¢. Thus in what follows, in the series representa- 


tion of the functions involved only terms up to and including ¢ will be retained. 
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In addition we also note that since (VX — 6/2) = —K/hy = —n'(K,/m) then with the 
restriction K/h, < u’, it follows from (2.7) that 
k=] 
and 
3. 
k ~— >? = Ss 
2 m 


By imposing the above stipulations, we now proceed to reduce the solutions given in 
Sees. 3 and 4 to those valid for shallow shells of revolution. 
Case A: 


1) b, 1, 0. 


From (3.2), 


und hence by the first of (2.8), (2.17) and (3.1) 


, . lp, suds a a 
W=¢' (1 — i. )(a + iy). (5.1) 
Also by (2.14), (3.5) and (3.6) 
w=14+%% 
Po 
and hence (3.8) becomes 
: 1/2 1 p, ' : ale 
W=¢6 1 — 7p @ |[AJ2,3(n) + BJ_2/5(n)], (5.2) 
‘ 0 
— Do Do |” 5 pi 5/2 - 
7=t sh m 1+=*-¢|¢”. (5.3) 
o do ho ‘ Po 
2) b, 0, f a qi Po 
From (3.9) and (3.2 
} } nen 
~ =f + 16, (5.4) 
ho ho 


where 


hy _ Be ( “:) - 
hoa y+ h.)” (5.5) 


Again by the first of (2.8), (2.17), (3.1) and (5.4), 


W = [ + ; (P: +3 i - 2) |e + iy) (5.6) 


and by (2.14), (3.10) and (3.12) 


1 Po h, Pi) 
S o- ion ; i eieet fi} 
ial (?: * ho 
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with the aid of which (3.13) yields 


; ' l 0 h, at 
Yas E 4 (2 oo PN |iay, ia) + Bl-.A0), (5.7) 
10 \qo ho Do/ 
.2 Do Do} 3 ito , h DP 2 . 
f= =1m | — + — 2 pip. (9.3) 
3 Jo ho 10 \qo ho Po 
Case B: b b, b ¥ —1. 


Because of the several possibilities for the form of the thickness variation discussed 
at length in the previous section, the determination of the specific coefficients for the 
series expansion of any one of these possible forms of h (which can be deduced in a straight- 
forward manner) will not be considered here. 

It follows from the first of (2.8) that 


W € ) ¢ [1 ‘a (@ ~ 6s [a + iy) (5.9) 
ho 2 do Po 


and by (2.14), (3.6) and (4.11), 


5.10) 


va 
S 
% 
ee 
=| > 
ing eel 
ite. 
+ 
nN | — 
an, 
bo 
srs 
Q SO 
ee 
— 


Hence (4.19a) becomes 


W=¢ ( ) [ - (2? =. )o |s NAS. code + B¥ caccacdel, Gan 
\ho f Po Jo 


where 7 is defined by (4.19b) and (5.10). 
When h is uniform (in which case, as discussed earlier, / is approximated by a con- 


stant), with the aid of (4.19b) and (5.10), (5.11) reduces to 


—— 1/4(2p,/pPo — 41/ Qo) 1 pV ~ 416 
W=¢6 {1 - (b/2 42) P fl Azasvsorn(0) + BY20+y/0+n(m)], 6.12) 


where 


_ ws Do (1 + db) |”? s2ys2 L (2B — su) 2+ oe 
n = | am h re, | re) [ +5 “ rs b/24+2° ; (5.13) 


As an example of the foregoing solution (5.12), we now consider the case of shallow 
paraboloidal shells of revolution, the middle surface of which is specified by 


2 (4), 5.14) 
\do 


where n(n > 1) is a rational number. Also, the principal radii of curvature of (5.14) are 


a n—-1)—*-1 ~ 46 
,, = see’ o(tan ¢) ; (5.15) 
n— 1 

a (tan @) 
n—2 sin @ 
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where a, as a representative length, is given by 


1/(n—1) 
Ao 
a= a . 
n 


Expansion of (5.15) in the series form of (2.17) leads to 


l 
b = ——], 
n— 1 
a : , =0 etc 
se ee " (5.16) 
Go = a, qa = 0 etc., 


the substitution of which in (5.12) and (5.13) results in the following solution for uniform 
shallow paraboloidal shells 


W = ¢'*[AJ2,,.(n) + BY:2,,(n)], (5.17) 


1/2 
n=? [am ‘i 1] eer, (5.18) 
ly it 


The above example of shallow paraboloidal shells of uniform thickness has been 
treated earlier by E. Reissner [6] and it would be of interest to show the correspondence 
of the two solutions. Using his own notation (except for a subscript R which, when 
necessary, will be added in order to distinguish the symbols from those of the present 
paper), Reissner in [6] defines the middle surface of paraboloidal shells by 


where 


r= ak, z= durf r(é), fr => é", (5.19) 
‘where a is a reference length, {; is of order unity and uz a number small compared with 
unity. A shallow shell is defined by the stipulation that for it, terms of order uz 
may be neglected compared with terms of order unity.”” With this stipulation, a = a, 
A —3/(4E°), uw? = m/ho a nu , Vi = &"-*, and (2.13) becomes (which, except for 


a constant, is the solution given in [6]) 


a 
We - £ os J An )( (5.20) 
1 Y.,.(n*)| 
where 
1/2 
g* =2 [am = as eet (5.21) 


Prior to establishing the relationship between the two solutions (5.17) and (5.20), 
we note that if terms of 0(¢°) are neglected in the series expressions of r and z of this 


paper, these quantities read: 


r=ag’"”, z=-¢9""". (5.22) 
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Thus, the form (5.19) is equivalent to (5.22) of the present paper if 
l ~ 6 
Mr (3.23) 
7 
‘ 1 
é QP 
It follows from (5.23) that. * n, and hence 
Vv é Wr (5.24) 
It may be noted that though the two solutions W and W, differ, due to the particular 


choice of the independent variable, they yield identical expressions for B + iy defined 


in (5.9) 
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INTERMODULATION PRODUCTS FOR y-LAW BIASED WAVE RECTIFIER 
FOR MULTIPLE FREQUENCY INPUT* 


BY 
EK. FEUERSTEIN 


Laboratory for Electronics, Inc., Boston, Mass. 


Abstract. The intermodulation products obtained by passing the sum of N + 1 
sinusoids of amplitudes P, , --- Py through a power rectifier of characteristic 
] \0 V<B 
laV—B) V>B »>0 


have been expressed by 8. O. Rice and W. R. Bennett in terms of contour integrals 
involving products of Bessel functions. In this paper these integrals are rewritten as 
improper integrals on the real line plus constant terms. These integrals converge fast 
enough in many cases to be useful in numerical integration. 

\ number of sub-cases arise, for which formulas for the intermodulation products 
are listed. 

Introduction. The evaluation of the intermodulation products obtained by dis- 


torting the time function 
N 
V(t) = pm P, cos (p,t + y,) (1) 
rough il nonlinear device 
I = af(V) (2) 
heen considered by several authors-in recent years [1-10]. 


Che output can be described ast: 


I(t) > ney > 7 ae ae a (3) 


mo=O0 my=0 


COS Mo(Pol + Yo) Cos m,(pit + ¥:) +++ COS my(pyt + yn), 


E, l, 
ée 2, m = 1,23,9-**. 
\ typical term of (3), expressed in terms of sum and difference frequencies, is 
Lngss+my COS [Mo(Pot + Yo) %& m (pit +1) & -++ My(pyl + vy]. 


It has been shownt that for a wide class of nonlinear functions the intermodulation 


oefficients are obtainable as contour integrals 


*Received April 25, 1956. 
Reference [3], Eqs. (4.9-15). 


tReference [3], Eqs. (4.9-17). 
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ad . N 
ee | F (iu) I] J, (Pw du, (4) 
T Je r=0 
where 
M= >) m (5) 
and 
F(iu) = | I(V) exp (—iuV) dV, (6a) 
KV) == | F(iu) exp (iV) du. (6b) 
2r Jc 


For the biased v-power rectifier for which 


= \0 V<B (7) 
la(V—B)’ V>B 
v> 0, 
F(iu) = al(v + 1)(iuy*” exp (—7uB) (8) 
and the path of integration C is the real wu axis from —@ to +© with a downward 


indentation at u = 0. 
Thus for »-law rectifiers, (4) specializes to: 


*M . 


N 
i , 
ia al'(vy + 1) | (iu) ‘’*"” exp (—iuB) [] J, (Pu) du. (9) 
7 v~f r 0 


In this paper the contour integral (9) is transformed into real converging integrals 
with bounded integrands plus a constant. Since most of these integrals converge fairly 
rapidly, numerical desk or machine computation of the coefficients becomes quite 
practicable, especially when only moderate accuracy is required. 

Notation. It should also be noted that the coefficients are functions of v, my 
m,--: my and P, , P - Py , the bias B, and the seale factor a, all together 2N + 5 


degrees of freedom. However, (9) is homogeneous, and can be written 


*M 


Ps, _(B\ Fr 
Pe ee aPsT(v + 1) . | (iu) °"*”’ exp ~iu(*) I] J AP uas/P) du. (10) 
T v¥¢ 0 ) 


It is convenient to choose as P, the largest amplitude of the set of input amplitudes, 


and to substitute 


— =k,<1, k=1, (11a) 
at an) 


The integrand of (10) becomes a function of v, the degree of the rectifier law; h, the 
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normalized bias; k, «++ ky the normalized input amplitudes; and my --- my the order 
of the harmonies; altogether 2N + 3 variables. 
Kquation (3) may be rewritten as 


x N 
I(t) = aPi dS +e SY FA ay(hy ky, +++ kw) T] en, cos m(p.t + y,), (12) 
where the coefficients 
a” Ae” cath, by ~++ Rg) @ AS? craze (13) 


The coefficients on the left [with the variables (h, k, --- ky) spelled out] correspond to 
the notation of Sternberg, Kaufman, Shipman and Thurston, who have called the 
coefficients “Bennett functions,” while the coefficients on the right are given in the 
notation of Rice, provided that the amplitudes have been normalized. This paper will 


use the notation of Sternberg, Kaufman, Shipman and Thurston, and concern itself 


with the transformation of 


-M 


. V 
1 (oe o« Kw) (vy + 1)! | (71) **) exp (—iuh) I] J...(k u) du, (14) 


“Cc 


M=) m,, k=1, 


into real integrals. The reader is advised that for integer v, especially of low order, the 
coefficients (14) are available in terms of published tables [6, 8, 9] computed by numerical 
evaluation of certain multiple integrals, and may be extended to other values of the 
parameters by recurrence relations [2, 4, 5]. The utility of this paper rests partly on the 
independent evaluation of coefficients without using multiple integration, partly on the 
evaluation of coefficients for non-integer power laws for which the multiple integral 
formulation does not apply directly, and partly on the function theoretical perspective 


concerning Bennett functions, obtainable from the integral formulation. 


Evaluation of the contour integral. ‘The integral (14) is broken up into three ranges 
of integration: over the negative real u axis, © < u < —46; over the positive wu axis, 
6 <u < ;and asemicircular contour with radius 6 indented into the negative imagin- 
ary half plane. By a change of variable from u to —u the negative w axis integral is 
reflected onto the positive half axis. One obtains: 


a 


" 
Ae svcmul hs By 9 *** By) = . Tv + ned f (iu) “°*” exp (—-7ul) I] J,, (kw) du 
T Le 0 


4 r 


~ «x 


N 
a aty42) a 
+ 2 cos 2 (vy + 1) | u~’’*? eos uh I] J, (ku) du (15a) 


3 


7) 


V \ 
—(v+1 . 
u~"*” sin uh | J (ku) du 


j 


° T 
— 2sin 5 + 1) | 


é 


? 
Zz m, = M_ even 
r=0 
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or 
| ; ( . N 
’ . -M . . y+ , 
De ae Ohl 5. 8° * ey (py + 1)0" "$2 p (7u) exp (—tuh) I] Jn (kw du 
T J 3 r=0 
»2 NV 
+ 2 cos 5 (vy + 1) | u~"’*? sin uh I] Jn (kw) du (15b) 
~ e r 0 
fa N ) 
+ 2sin a V + 1) u-""” eos uh I] Jim (ku) du ( 
if 


>. m M odd. 
1 by 


The semicircular contour integral in (15) is evaluated by a) replacing the integrand 
a power series in u, >) replacing u by its complex polar form 


u 6 exp (20 du = tu dé (16) 
and c) integrating term by term from @ r to 6 2m. 
The resulting series in ascending powers of 6 will contain negative powers of 6 (the 


principal part) provided pv M; a series in ascending positive powers of 6, and, if v is 
“a, , Where 


an integer, a term free of 6 (the zero power of 6) equal to m 


. [exp (—iuh) [] J,,,(k,u)] (17) 


Vv. Ldu 


the residue of the integrand of (14) and is obtained for com- 


a, can be recognized as 
putational purposes most easily by identification with the vth coefficient of one of the 


power serie ~ expansions 


Fy(u cos uh I] J (kw) > a,u'; if M-+vp odd (18a) 


N 


Fu isinuh [] Jn(kuw) = Doau'; if M+v_ even. ISb) 


Since all powers of u have zero coefficients for 


ee M in the ease of (18a) and 


0<I1< M+41 inthe case of (18b), 


the residue a, vanishes for M > v + 1if M + vis odd and for M > vif M + vis even. 
The real integrals of (15), containing 6 as the lower limit, can be transformed so as 

to bring into evidence the principal part of the Laurent series in 6 contained in them. 

If F(u) and its derivatives F’(u), F’’(u), F‘”’(u) are bounded and p is an integer 

so that 

p-—l<v<p (19) 
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then, by repeated application of integration by parts, one obtains the identity 





i Fw) ay i [xo ,—_F) ee ‘fie, =a | 
5 ou vb viv — 1)b vy—1)---v’—pt+)b” 
Ws Ald \ AAS 1) 
a [ Fe) 4 . ee per) | (20) 
v6 viv — 1)d vy — 1)---@~—-pt+)ar” 


b 


—(y=pt1) i 
ig vv — 1) ---@v’—pt]) i u (u) du 


Substituting for F(u) F,(u) or F,(u) (see 18), one sees that the first bracket of (20) 
vanishes as b — . Replacing each term of the second bracket by its power series ex- 
pansion, one obtains a Laurent series in 6 for the second bracket collectively. 

From (15a) and (15b) it is clear that for v integer wu ‘’*'’F(u) in (20) must be even 
for which u “F(u) is an odd function. Hence, for all cases of interest in which » is an 
integer, each term of the brackets of (20) will be an odd function and it will be free of 
the zero power of 6. If v is not an integer, all powers of the Laurent series in 6 of the 
second bracket of (20) are not integers; therefore a possible term constant in 6 is excluded. 

Finally it can be shown that the integrals on the right hand side of (20) contain, as 
far as terms of interest for substitution into (15) are concerned, 6 only as a power series 
in positive terms of 6. Thus, as 6 becomes small, these real integrals converge to their 
value with the lower limit equal to zero. 

In the case of integer v, » = p and the integral on the right of (20) becomes: 


1 ss) 


vi Js 


uF’ (u) du. 


From (18) u’F’’~'’(w) is an odd function; it follows that u~'F'‘’’(u) is an even function, 
and its integrand becomes an odd function of the limits of integration. On the other 
hand, since F''"’(u) exists for u = 0, its power series expansion involves only non-negative 
odd terms, the lowest of which is u'. Hence, the power series expansion of u~'F'‘”’(u) 
contains only non-negative integral powers of u. Integrating term by term and intro- 
ducing the lower limit one obtains only odd integer positive powers of 6. 

If vy is non-integral, the integrals on the right side of (20) involve the lower limit 
5 only in positive non-integral powers of 6. Since 0 < » —p + 1 < 1 and F’? (u) is 
expressible in a power series in non-negative powers of (wu), a term-by-term integration 
leads to positive powers of the lower limit only. We note in passing that, for the special 
case of M 0, 


Sa) 


Lim [ u"?*? eos uhd (uw) J (kyu) «++ Jolkyu) du 

weg A 
is an improper integral which for machine purposes may be transformed into a proper 
form with bounded integrand by one further integration by parts. 

Since the integrand of (14) is regular for all finite u, except at the origin, (14) must 
be independent of the value of 6. This implies that the principal parts of the Laurent 
series in 6 of the expressions in (15) cancel exactly. The contribution of the positive 
powers of 6 in each expression of (15) can be made as small as desired by making 6 small. 
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If (20) has been substituted into (15), the integrals on the right side of (20) go over into 
integrals with lower limit of integration zero. 


The only possible contribution of the semicircular contour is the residue term appear- 
ing for integer v 


The integrals (15) may be written: 


| Tt a f ad’ > 
| CO 9 (yp { | : | u — du | cos uh IT . ae ca | du 


21a) 
ie |2 5 T Ae : 
| in ; (vy 4 » | : | eile oe sin uh I I Jak) | du 
pit a, > if M even, 
, ; ‘(pv 1) Ww! 
{ } - | l 
vy l) +++ G p + 1) 
af 1’ . 
| cos ~ (p | » | | eo = sin uh [] 1.820) | du 
2 wT. ad’u oe 
2b) 
| - 9 . ‘la V 
| sin 9 (vy + 1) c | et a cos uh I [ J mA(k,1) du 
p! 1 “a, ° if V odd 
and for both Y Ila and 2 1b) 
0 vy non-integer 
a ) \ (21c) 
; exp tuh | | J nm (k,u) vy integer 
du ; 


as well as: 


Ko Lg i ie Es 
Equation (21) is more complicated than necessary since all possible cases are accounted 


for in two possible forms. Specializing further, one obtains simpler equations which are 
listed in Table 1. 


If 
N 
h> > k.: Oe vicgstily Me: ky) = 0 (22a) 
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N 
h<— Dk,: Ae? .comy(hy Rr °°* Ky) = Qe o*~ a, (22b) 


and y integer where a, is given by (17) or (18). 

The above results were obtained by observing that if v is an integer one may, without 
contributing to the integral (14), close the contour over the negative imaginary halfplane 
if h > y ey » k, ; and close the contour over the positive imaginary halfplane if h < — 
>>”, k, . These forms express the output when the bias is either so large as to prevent 
all output, or so negative as.to be free of cut-off regions. 

Concluding remarks. This paper has been concerned with the analytic formulation 
in terms of integrals with infinite upper limit. A glance at the original integral (10) 
should convince the reader that the integrals resulting from it must be bounded for all 
portions of the path C except possibly in the neighborhood of the origin, and it has been 
the burden of this paper to show that the contribution of the integral (10) over a portion 
of the path near the origin is also bounded. 

For numerical evaluation of the integrals in Table 1 one must either find closed 
form solutions in terms of other transcendental functions for which tabulated results 
exist, or resort to numerical integration. For the zero bias case with two frequencies one 
can indeed identify integrals (21a, b) with Weber-Schafheitlin integrals [1, 11, 12] provided 
that m) + m, > + 1. These integrals are expressible in hypergeometric series. For 
other two-frequency coefficients and integer power laws, the evaluation of the inter- 
modulation coefficients might best be performed in terms of the published Tables of 
Sennett functions, and recursion relations listed in [4, 8]. (These relations and others 
ean, of course, also be obtained by application of the standard recursion formulas of 
Bessel functions to the integrals in Table 1). 

In all other cases the integrals (22) in Table 1 must be evaluated numerically after 
truncation at an upper value of u. From the point of view of precision computation it 
is imperative that the incomplete (truncated) integrals converge to their limit rapidly 
as the truncation value is increased. Therefore, for large u, the power of u in the integrands 
(22) should be as negative as possible. 

For integrals (22c to 22h) the identity (20) can be applied in the form suggested in 
note 4 of Table 1. The integrals (22) are split into two regions at u = 1 and the slowly 
converging integral in the region of large u is replaced, by virtue of (20), by an integral 
with faster rate of convergence plus an expression in terms of the derivatives of the 
integrand evaluated at u 1. (See note 4 of Table 1). The problem of truncation still 
remains. Dr. Sternberg has pointed out to the author that, for | x | > 1 and | x | > m, 
the Bessel functions may be replaced by one of their asymptotic representations, (11, 


13). The simplest representation becomes 
. _ 1 . 1/2 ~ . i= Lud 9 
J (x) = ($42) cos | x + (m — 5) 5 ; (23) 


The above approximation can be applied to the estimation of the integrals (22) 
beyond their truncation value. Using x, = k,u, one obtains integrals involving trigono- 
metric funetions and negative powers of u. If it happens that the negative power of u 
is an integer, one can, by repeated integration by parts, obtain integrals with integrands 
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u~* sin u or u' cos u. Thus estimation of error resulting from truncation at U amounts 
to table look-up of Si(U’) and Ci(U). In general, an upper bound to the approximation 
is found by neglecting the oscillating cosine term. Unfortunately this bound is too large 
to be of much use in accurate computation. In general, the error in the truncation 
estimate, which results from a substitution of the trigonometric functions for Bessel 
functions, can also be estimated from consideration of the error in the first neglected 
term of the semi-convergent series (23). It should also be noted that, for small /, , u 
must be large to make k,u > 1. The latter fact precludes the direct use of the asymptotic 
approximation for the Bessel function in those cases where the input amplitudes P, 
[Eq. (1)] differ greatly; however, in this case one may use the power series expansion of 


the Bessel function for the intermediate region: 
U <u <- , A> 1. 


As for the numerical integration of (22) between fixed limits, the reader is referred to 


books on numerical analysis for a discussion of methods and error bounds. 
Acknowledgement. ‘The author wishes to thank Dr. R. L. Sternberg and Mr. J. 8. 
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— NOTES — 


LIMITS ON THE ZEROS OF A NETWORK DETERMINANT* 


By A. PAPOULIS (Polytechnic Institute of Brooklyn, Brooklyn, New York) 


The properties of a network are determined from the location of the zeros of the 
network matrix and its minors; the exact determination of these zeros often involves 
the solution of a high order equation. In many applications, e.g. in questions of stability, 
one is interested not in the exact determination of these zeros but in their exclusion 
from certain parts of the plane. In this note, limits will be given on the location of the 
roots of the matrix of a network consisting of two kinds of elements. 

The matrix of a network containing active sources and only two kinds of elements 
is given by : 

}ay.z + bir, Gye + Oe, tee 5) Aine + bin| 


aim [Ot i Sethe os re, r 
LdmiZ + Dat» Gn2® + Dany ey Anne + Dard 
where 
a.>0, bs>QO. (2) 
In this matrix, z = p for the R-L and R-C case, z = p’ for the L-C case. 


We introduce the constants’: 


R, = Zz Tie ; P, = > Sg b,, 


— a ; 
8, = Aresin —, yg, = Aresin— , (3) 
ay; Dex 
and 
{ Dy + P, in by — P, 
; max. — — a= man. ———, 

l Qk: R, 2 Qu tT R, 

a = max. (0, + ¢,). (4) 
" 


Suppose that z; is a root of the determinant of 1/(z). We shall show that z,; satisfies 


the following inequalities 
(5) 
Consider the system 


pe (a,,2; + b,,)w, = 0 oe Pee (6) 


*Received April 6, 1956; revised manuscript received June 5, 1956. 
‘In this paper we shall denote by 5°} ax, the summation for r = 1, 2, -+: ,nr ¥ k. 
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The determinant of its coefficients is by assumption zero, therefore it has a non-trivial 


solution w, ,w.,°+** , w, . If w, is the variable with the greatest absolute value, then with 


we have 
Ee F r=1,2,--- ,n. (7) 


Dividing the kth equation of (6) by w, and solving for z; we obtain 
| “ 


bu + Do’ dart, 


-2, > owe’ (8 
tat 2 Gute 
Therefore [see (2) and (7)] 
ba Pe ey, cbt Pe 
Qu. + R, ~ a, — R, 
| arg (—z,) | < Oe + ¢- (9) 


The relations (9) are true for some k, hence (5) follows. Interchanging rows and columns 
of M(z) we obtain results similar to (5). Combining the two sets of inequalities for z, 
and arg (—z,;) we might thus improve our limits. 

From the above it follows that if 


T 
(10) 


then 


oe ae ai 
) aE™(—4%i) | Qo 


hence z; has a negative real part. This is true if 


Axi b 
a = De kk _ ae eens (11) 
R (Q)'/? ; is ™ (2)'/ k Lb > . LI 
because then 
a wt “ T 
oS ’ Pi . 
{ 4 


These conditions can readily be verified for a given matrix; one thus establishes simple 
sufficient conditions for the negativeness of the real part of the zeros of M(z), which is 
useful information in network analysis. Suppose that a network contains R-C elements, 
control sources and feedback, depending on a parameter k. The matrix of such a network 
is given by (1) and the system is stable if the real part of its zeros is negative. To determine 
the values of k for which this is true, it would be necessary to find the zeros of the de- 
terminant in terms of k; this is not easy. The above discussion offers a simple method of 
establishing sufficient conditions on k, for the stability of the system; clearly beyond 


these limits the system is not necessarily unstable. 
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HEAT TRANSFER AT RECTANGULAR CORNERS* 


By J. C. WILSON (Research Laboratory, Redstone Arsenal, Huntsville, Alabama) 


Introduction. In this paper, a solution to the heat equation, Au = u, , appropriate 
to a small neighborhood of a rectangular corner is developed. It will be assumed through- 
out that the initial condition is that of zero temperature. The boundary values will be 
given by a continuous function of position. The equation is solved by means of “Green’s 
function” and the solution expressed as an infinite series. This series is relatively simple 
to evaluate and the error involved is readily calculated. Moreover, the solution may be 
differentiated with respect to any one of the variables involved. 

I. Let the rectangle be such that a vertex is at the origin of a (£, 7) coordinate system 
with its boundary consisting of the positive — and 7 axes, together with two half lines at 


infinity. The boundary function is defined as 


kill — €/a’] 7 = 0, O<s£ < a, 
v(f, nn) = kf — n a’}, é == 0, 0 4 n < a, (1.1) 


0, for all other points (&, 7) on the boundary, 


where a and k are positive, real numbers. From a physical viewpoint, the function ¢ 
may represent a good approximation to boundary functions resulting from various 
modes of heating, at least for — and 7 small. For example, suppose the source is a two- 
dimensional flame of arbitrary outline, applied to the corner in such a way that its 
temperature distribution is closely approximated by either T, = k [1 — (€ + 9°)/a? + 
((én)], where f is a polynomial in &y, or T, = k [1 — (€ + 9”)/a° + (higher powers of 
£ and »)}, with 7, and 7, being positive or zero for all £ and 7. If the flame temperatures 
are best approximated by 7, , then for0 < — < a,0 < n < a, the boundary temperature 
must be equally close to ¢. If T, is best, then, for & and » small enough, one has boundary 
values again arbitrarily close to ¢. 

II. A solution of Au = u, will now be obtained which is valid for x, y, and ¢ positive 
and such that 

(a) limuw=0, for x>0, y> 0; 


(b) u = ¢ fort > 0 and (x, y) on the boundary. 
First a “‘Green’s function” is defined as follows: 
G(x, y, t:&, 9, 7) = A(t — a) fexp [—( — a)’/4(t — 1] — exp [—€ 4+ 2)’/4(t — D)]} 
-fexp [—(n — y)’/4(t — 7)] — exp [—(n + y)*/4(t — D]}. (2.1) 


{s is well known, by employing Green’s formula, one arrives at the solution 


u(x, yy, = / | g(a, y) = as| dr, (2.2) 
0 7R 


’ 


*Received April 16, 1956; revised manuscript received June 26, 1956. 
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where the symbol {, denotes the integral taken over the boundary of the rectangle 
and 0G/dn; signifies differentiation in the direction of the inner normal. After sub- 
stitution, differentiation, and simplification, (2.2) becomes 

- 


u(x, y, t) = y/m | g(é, 0)[exp {(—20[E — x)? + ¥/[E — 2)? +’) 


aa 


— exp {(-i0)(E +a’? + YN /E+a° 4+ J) dt + 2/r | yg(0, n) 


| “0 (2.3) 
[exp {(—30[(n — y)? + 2°)}}/[(a — y)* + 2°] 
— exp {((-2)[(n ty +r }/aaty’ +e )yldn=L-Lh+h—-h, 
where a is the same as given in (1.1). In order to show that lim,... w = 0, consider the 


integral J, . The integrand may be expanded in a series converging absolutely and 
uniformly for all x, y, and t bounded away from zero. Integration yields the new series, 


aa 


be y/n exp (—y"/4t)-1/n! 1/(40” | g(€, OE — x)"/[E — x)” + yy’) a, (2.4) 


n= J0 


convergent under the same conditions. Now the nth partial sum of (2.4) has limit zero 
when ¢t — 0. Thus, the limit function approaches zero with ¢ so that J, — 0 as t — 0. 
The same result holds for J, , J; , and J, . 


To show that wu takes on the prescribed boundary values, consider the expression 





dI, = y/n exp {(—40[€ — 2)° + y'}S/[E — 2) + y'). 


For x + &, lim,... dJ, = O and for x = £,dI, — ~ asy — 0. Also, lim,... {2.. d/, dé l, 
so that dJ, has the “character of a 6” function. Using Schwarz’ distribution theory, 


together with the fact that ¢(é, 0) can be extended to be C™ for all &, the limit 


vr 


lim | g(é, 0) dJ, dé = g(x, 0). 


yo0 Ja 


But for x > 0 and bounded away from zero, 


lim | g(€, 0) dl, dg = lim | (€, 0) dl, dé 
and this limit is uniformly taken in the sense of the Moore-Osgood theorem. Also, for 
y > 0, 
lim | g(é, 0) dl, dé 


*zro YO 


clearly exists for all > 6 > 0. Therefore, 
lim | g(t, 0) dI, dé = g(x , 0), 


for all x) bounded away from zero. It follows immediately that 


lim [-—/J, + /,; — 4,] = 0. 
Cs.) 
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Using identical arguments, 


lim | 


0 /0 


Gas 


é g(0, n) dl; dn = (0, yo) forall y > 6>0. 
Again, the remaining double limits are zero. It remains to show that 
lim u = (0, 0). 
(J=0 ) 
For this purpose, it is sufficient to consider the sum 
R-h+ii3-l%, (2.5) 
where J, isJ, ,a@ = 1, --- , 4, with the exponential term set equal to 1. After suitable 


transformations, one can write (2.5) as the sum 


2/mig(x + yA, 0) arc tan x/y + 9(0, y + 2X’) arc tan y/z} 


Lf f°?” oa + yo, 0) “9 oa — yo, 0) 4 | 
= J, ie I Ito “YS 
af pee 0, y + x0) f **»/* 00, y — xa) S om 
Fo l+oa ee ae l+o ” eine aca: 


where —x/y <A < Oand —y/x < ’ < 0. Since ¢ is a continuous function of z and y 
in any neighborhood of the origin, it is certainly possible to find a 6 > 0 such that for 
r < 6, y < 6, the following inequalities hold: 

A — ¢(,0)|<«,|B|<e, and |C|<e, for ¢€>0O. 
The existence of the limit in question is therefore assured. 

III. To obtain the series development of the solution, first the exponentials in each 
integrand of (2.3) are expanded. The resulting expressions are series converging absolutely 
and uniformly for x, y, and ¢ in any closed, bounded domain, minus the origin. After 
substituting for ¢ as given by (1.1) and collecting similar terms, one may integrate, 
obtaining the expression 
( k/ra?{1 — (a° + y’ — 2°) are tan 2ry/(a’ + y’ — 2’) 

a +2” — y’) arctan 2ry/(a’ + 2” — y’) 

ry In (a? + y’)*/[a® — 2a*(x* — Bxr*y? + y*) + (2? + y’)‘)} , 

(3.1) 
— ka’xy/16rt? + ka’xy/96rt*[a’?/3 + 2” + y’] 


9 


ka*xy/3072rt'[a*/2 + 8a°/3(x? + y*) + 3(2? + y’)’] 
+ ka°xy/30,720rt’[a°/5 + 5a‘/3(a2? + y’?) + 10a7/3(2? + y’)? + Ae? + y)*] -— --- 


where the principal value is used in computing the arctan. 

The series is alternating and converges absolutely for all x, y, and ¢ positive. For 
applied purposes, one may have to restrict the values of f to some definite range in order 
that the solution give results in agreement with actual conditions. Appropriate ranges 
for ¢ would, of course, be governed by the nature of the heat source and the material 
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heated. In any case, it is clear that the most rapid convergence is obtained for small 
x and y and comparatively large ¢. Although the first term appears formidable, it turns 
out that in many instances the arctan and log terms are negligible. For example, if 
k = 10°, a = 2, and zy < 4/10", the arctan and log may be neglected if an error within 
+1/10 is satisfactory. Since the series is alternating, the error involved is less, in absolute 
value, than the first term not employed in the calculation. Thus, once k and a are fixed, 
one might check terms, beginning with the third, ie., ka°xry/96mt* [a’/3 + x7 + y'], in 
order to obtain a range for x, y, and ¢ so that a given error tolerance is not exceeded. 
From the form of each term, an inequality involving the product, xy, is simplest to 
handle and furnishes a good check for x and y small. As a usual occurrence the calculation 
of three or four terms gives sufficient accuracy for applied purposes provided a reasonable 


balance is maintained among the variables 


VISCO-ELASTIC STRESS ANALYSIS* 
By J. R. M. RADOK (Brown University) 


1. Introduction. In his paper on stress analysis in visco-elastic bodies [1]** E. H. 
Lee bases his reasoning on the concept of an associated elastic problem to which a 
visco-elastic problem reduces after removal of its time dependence by application of 
the Laplace transform. Thus Lee’s method requires the application of the Laplace 
transform to the boundary conditions as well as to the basic equations and it might be 
expected that it is restricted to problems whose boundary conditions admit such an 
operation. As a result, for example, the problem of indentation of a half-space by a 
curved punch could not be solved by this method, since at any one point of the boundary 
at different times stresses or displacements are specified. 

It is the purpose of this paper to extend the applicability of Lee’s method to problems 
of the above type and to show that the apparent restriction is due to the process by 
which Lee deduced his method, in particular, due to the concept of the associated elastic 
problem. 

At the same time, the Laplace transform method will be restated in a different form 
which may be called the method of functional equations. This method is completely 
equivalent to Lee’s method, since both these methods coincide, if the functional equations 
are solved by operational methods. However, the extension of the applicability of the 
Laplace transform method to the wider range of problems requires the functional equa- 
tion approach for its justification. 

2. The method of functional equations. The basic, quasi-static equations governing 
the linear theories of isotropic, elastic or visco-elastic media, referred to orthogonal, 


rectilinear coordinates x, , may be written in the form 


1 (du; ou 
Tet qa r Ys i i 
P’s:, = @'e;; , R'o;; = Se; é; = lilies ae (2.2) 
ei -s ; 2 \dz; dx,/’ 
*Received April 25, 1956. The results presented in this paper were obtained in the course of research 
sponsored by the Office of Naval Research under Contract Nonr-562(10). 
**Numbers in square brackets refer to the bibliography at the end of the paper. 


1957] J. R. M. RADOK 199 


where o,; , € ;; , u; and X, are the stress, strain, displacement and body force components, 


s;; and e;; are the stress and strain deviators, defined by 


] 


83; = ¢y - 46445;; ’ C3; = i. — 46,45; ; ’ (2.3) 
and p’, Q* , R’ , S* are linear differential time operators with constant. coefficients 

Pp 3” @ 3” 

P=Yasm, V=Lap 

2 Pn 34 Y= de We 5 


s Ee, e 
S = Lis op 


In the particular case of an elastic material, the orders p, g, r, s of the operators in (2.2) 


Rr’ = p .. < . 
(2.4) 


are all equal to zero and hence 


So 


fo Mm, = 3k, (2.5) 
Po To 
where » and k are the shear and bulk moduli repectively. 
Any problems in the theories of elasticity and visco-elasticity will be completely 
specified, if the above equations are accompanied by boundary conditions which mostly 
take either the form of given tractions T’, 


o;n; = T; on L (2.6) 


2 
or of given displacements f, 


u; =f; on L (2.7) 


or of a combination of these conditions, where L is the boundary of the body and n; 
are the direction cosines of the outward normal to L. In addition, in general, due to the 
time derivatives in the stress-strain relations, initial conditions have to be satisfied 
the number of which depends on the orders of the operators P’,Q*, R’, S*. In the case 
of an elastic material no initial conditions occur, since p = q = r = s = 0. 

There are different ways in which time dependence may enter into the solutions of 
the basic equations (2.1), (2.2). The given tractions 7’; and displacements f; may be 
functions of the time, the sections of the boundary to which the conditions (2.6) or (2.7) 
apply may vary in time or a combination of these conditions may occur. In general, 
the solutions will thus depend on the time ¢ as well as on the space coordinates x, and the 
material constants p, , Qn, tn, Sn - 

First of all, it will now be shown that every elastic solution of a given boundary 
value problem, for given initial conditions, uniquely determines a visco-elastic solution. 
The proof of this statement is effectively contained in Lee’s paper [1]. In fact, Lee uses 
this result in his second method when dealing with moving loads. 

Assuming for simplicity that the visco-elastic solution is subject to zero initial 
conditions and that therefore the given surface tractions and displacements, specifying 
the boundary value problem, vanish for ¢ < 0, the time dependence of the basic equations 
(2.1), (2.2) may be removed by application of the Laplace transform. Denoting by « 
star the Laplace transforms of the time dependent quantities in these equations, one 
finds in this way 

oot +. Xe = 0, (2.8) 


Ox; 
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P*’s* = Q*%*, R*'o* = S*'e* (2.9) 


where for a visco-elastic material the operators P’, Q‘, R’, S* become polynomials in the 


Laplace transform variable. Since (2.9) also refers to the elastic case when 


P® =p, Q*=qm, R* =n, S*® =H, (2.10) 


the Laplace transform of a visco-elastic solution may be obtained from that of the 
elastic solution by means of the substitutions (2.10), when not all of the orders p, 4, 
r, s of the operators are equal to zero. However, in many cases, the elastic solution will 
permit a rearrangement of the elastic coefficients po , Go , 70 , So , So that the Laplace 
transform of the visco-elastic solution is equivalent to ordinary differential equations 
in time with initial conditions. 

This result justifies the following alternative procedure for deducing a visco-elastic 
solution from a given elastic solution. Instead of applying the Laplace transform to the 
elastic solution, the elastic coefficients pp , % , > , 8 in this solution may be replaced 
by the operators P’, Q’, R’, S°. This process leads to functional equations for the stresses, 
displacements, etc. of the visco-elastic solution. In many cases which involve the elastic 
constants in a simple algebraic manner, these functional equations will be differential 
equations in time which may be integrated for given initial conditions. 

So far nothing has been said about the boundary conditions, satisfied by the visco- 
elastic solution, obtained in the above manner from the elastic solution. If the boundary 
conditions of the elastic solution can be transformed by means of the Laplace transform, 
i.e., if an associated elastic problem in the sense of Lee exists, both the elastic and visco- 
elastic solutions satisfy the same boundary conditions. 

However, it will now be shown that the above visco-elastic solution always satisfies 
the same boundary conditions as the elastic solution from which it was obtained. This 
may be done using the method of differential equations, irrespective of whether the 
associated elastic problem exists. For this purpose consider o,;n; or u,; or both of these 
expressions for the elastic solution at every point of the region occupied by the body, 
depending on which of these quantities enter into the boundary conditions. The first 
of these expressions is to be understood in the sense that inside L the direction cosines 
n; may be given any desired values. Since o;; and wu; , referring to the elastic problem, 
depend on the elastic constants as well as on the space coordinates x, and the time f/, 
the above quantities will depend on the coefficients py , Go , To , So everywhere inside the 
boundary L, except possibly at certain points or on preferred planes, such as planes of 
symmetry where shear stresses vanish. However, on the boundary L, where the values 
of the direction cosines are given, the above expressions, when they enter into the 
boundary conditions, do not contain elastic constants, since the tractions and dis- 
placements are independent of the properties of the material. Thus, the expressions 
o,;n; and u; lead by the method of functional equations to differential equations at 
points inside the region, occupied by the body, but not on the boundary L, wherever 
tractions or displacements are prescribed. Hence, the visco-elastic solution satisfies the 
same given boundary conditions as the elastic solution. 

It should be noted that the above result does not mean that the elastic and visco- 
elastic solutions agree all along the boundary, since, for example, at points where dis- 
placements are prescribed, there will be functional equations for the visco-elastic stresses 
which therefore, in general, will differ from those of the elastic solution. 
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In the next section, the method of functional equations will be illustrated by two 
examples; the first which deals with the point force on a half-space has also been con- 
sidered by Lee [1], since it has a corresponding associated elastic problem, the second, 
the two-dimensional problem of a varying size circular hole in an infinite body, subject 


to internal pressure, has no associated elastic problem. 


3. Applications. For the purpose of comparison with the Laplace transform method 
consider the case of a point load on a half-space z > 0. This problem was also used by 


Lee [1] as an illustration. The elastic solution is 


ae P(t) {( oe 2»| 1 = 5 (r? rs 2”) = — 3r°e(r? a eyw 


Qn 
By (2.5), 


9, — —_39o%o 
lata 2soPo + Too ’ 


and hence the corresponding functional equation is 


r 


QS'P” + R'Q)o(r, 2,1) = {sqrn'P(o| 4 _ 
— (28°P” + R'Q)P()3re(r* + at, 
For the Voigt material, considered by Lee, 
P=1, Q@=AS+B, R=1, S=C 


and (3.3) gives the differential equation 


; a | ee 0 1 Zope 2 aad 
( 20 +B+A 5/2 (r,z, t) = on {a( 4 at + B)PCo| + a (r° +2) 


—_ (2c + B + A 2 Pyare’ + ah. 


Solution of Eq. (3.5) for zero initial conditions and 
0 for t<0 
P(t) = 
P, for t>0 


by operational methods leads to the Laplace transform 


. P, f 3(Ap + B) 1 _ £72 2 =| —— Qy2 ofp? 2-5 + 
o* = Onn \2C + B+ Ap E ny +) sre + 2) 


Im 


= 


of the visco-elastic stress ¢, , and hence to the result given by Lee. 


(3.1) 


(3.2) 


(3.3) 


(3.4) 


(3.5) 


The problem of the moving load, likewise considered by Lee, is easily solved in the 
same manner from (3.5) after appropriate changes on the right hand side of that equa- 


tion. 


As another simple example consider the plane problem of a circular hole in an infinite 
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plane, subject to uniform pressure —P(t). The elastic solution gives for the radial 
displacement 
P(t)R® 
= —— ; (3.6) 
2Qur 
hence, by (2.5), the functional equation is 
. | dg pss 
Qv, = — P*P(b) 
. 

and for the same material as before, specified by (3.4), this equation becomes 


/ 4 \ >2 
(A 4B, = he P(t). 3.7) 
4 c rT 


Thus, the ge neral solution for v, is 


4 Bt sag . 
_ | e®/AD(t) dt + De®", 3.8) 
Ar 4 


ow 


where D is an arbitrary constant to be determined from the single initial condition 
required in this case. 
Now let also the radius R of the hole vary with time. The elastic solution (3.6) still 


applies and the boundary condition is now 
O, —Ii(t) for r = RU. 


This boundary condition cannot be made time independent by application of the Laplace 
transform. Nevertheless, by the results of Sec. 2, a visco-elastic solution is given by 
(3.8), where now R? must be included under the integral sign. 

4. Conclusions. It has been shown that the elastic and visco-elastic solutions 
of a given boundary value problem are related to each other through functional equations 
for the quantities specifying the solutions. If the elastic solution of the boundary value 
problem is known, these functional equations can be obtained by expressing the elastic 
constants in this solution in terms of the differential'time operators in the general 1so- 
tropic visco-elastic stress-strain law. In many cases these functional equations turn 
out to be ordinary differential equations in time whose solutions for given initial con- 
ditions make up a visco-elastic solution which satisfies the same boundary conditions 
as the elastic solution 

If these functional equations are solved by use of the Laplace transform, the later 
stages of obtaining the visco-elastic solution coincide with the Laplace transform method 
of Lee [1]. The deduction of Lee’s method was based on the concept of an associated 
elastic problem, requiring the application of the Laplace transform to the boundary 
conditions as well as to the basic equations. The present paper establishes the independ- 
ence of the Laplace transform method from the associated elastic problem and hence 


widens its range of applicability. 
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ON THE FINITE STRIP PROBLEM 
By CARL E. PEARSON (Harvard University) 


1. Introduction. Many two-dimensional problems of physical interest reduce to 
the problem of solving a singular integral equation, with finite limits, of the form 


, 


| [Pia —t)inix-—t!|+ Qa — Ojf(d) dt = g(a), e<e <6, (1) 


where P, Q, g are prescribed well-behaved functions. For example, let a perfect fluid 
move with constant velocity parallel to the z-axis. If the portion (a, b) of the x-axis is 
occupied by a thin metal strip heated to a prescribed temperature g(x), then the temper- 
ature distribution in the fluid is easily obtained if the heat-source intensity f(¢) in the 
interval (a, b) is known. This quantity f(¢) satisfies an integral equation of the form 
of Eq. (1) where the kernel (quantity in square brackets) is now proportional to 
lexp a(x — 1) Ko(a@| x — ¢})]; here Ky is the modified Hankel function and @ is a constant 
involving the properties of the fluid. 

No general closed-form expressions for the solution of Eq. (1) are known. In practice, 
three different techniques of approximation are used, depending on the length of the 
interval (a, 6). For intervals (a, b) so long that the strip may be considered semi-infinite, 
the solution f(¢) near each end of the strip may be obtained by the Wiener-Hopf method 
|] and these two end solutions may be superposed to give a practically useful solution 
for the entire strip. If on the other hand the interval is so short that only the zeroth 
order term of P(« — ¢t) need be considered, use can be made of Carleman’s formula [2] 
for the solution of 


sl 


| f(t) In| a — t | dt = g(x) (2) 
* 1 
vhich is 
f(a 2 7 3 ' g'(s)-Q — s)' : ds = _ [ Li) ~ 2° (3) 
wr (1 q ‘ EP s—2 w(l—2)* M2J_,(1 -—s) 


Here and in future integrals of the form of the first term of Eq. (3) are to be interpreted 
in the principal value sense. 

The in-between case of an interval of moderate length is always troublesome and 
usually requires either infinite series expansions or numerical techniques. In the following, 
a method wil! be described which gives an exact solution (in closed form involving simple 
integrals) of Eq. (1) for that case in which P and Q are nth order polynomials. The 
method involves the solution of a set of n linear equations* in n unknowns and so becomes 
tedious as » increases; nevertheless, its advantages of exactness and straightforwardness 
of calculation should make it useful in certain “‘in-between”’ finite strip problems. 

Before proceeding, some general comments concerning Eq. (1) may be of interest. 
lirstly, if all terms except ff(é) In | « — ¢| dt are placed on the right hand side and 


Received April 25, 1956. 
*The coefficients of these linear equations are expressible in terms of Bessel functions at worst, 
which are easily obtainable from tables; thus no numerical integration is necessary. 
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temporarily incorporated into g(x), then the use of Eq. (3) with this modified g(x) 
results in a non-singular Fredholm integral equation for f(x). This method (due to Carle- 
man [3]) is of theoretical value but the resulting Fredholm determinants are too compli- 
cated for practical use. Secondly, attempts to extend the Wiener-Hopf method to 
equations of the form (1) have often been made; although the various Fourier transforms 
involved can be determined within entire functions, the essential difficulty seems to lie 
in meeting the requirement that an inverse transform vanish over selected portions of 
the real axis—a condition which is automatically met in the Wiener-Hopf method but 
not here. Thirdly, an interesting approach suggested by Latta [4] utilizes the differential 
equation satisfied by the kernel function to obtain a differential equation for the solution 
in those cases where g(x) is a simple function (such as a constant or power of x). The 
differential equation may then be solved by series expansions. 

2. Solution of Eq. (1) for polynomial P and Q. For convenience, the limits of 
integration will be chosen as —1 and +1. Then define, for any complex point z not in 
the strip (—1, 1), the analytic function S(z) by 


1 


.@ ve f ; fz—t ; 1 
S(z) = (2° — 1)” | Piz — t)In Cx ‘) + Q(z — ) |S dt, (4) 


where (2? —1)! has the branch line (—1, 1) and behaves like +2 for large z, and where 
In i: _ ) In ( t) In (zg + 1 
= zZ- _ Zz ). 
2+ 


Here and in future the value of In (z + 2x), real x, is chosen as real for z on the real axis 
to the right of the point —2; the branch line runs from —z to —o. 

Clearly, S(z) is single-valued outside the strip. Also, the limiting values of S(z) 
exist as z approaches a point x on the interval (—1, 1) from above or below, and are 
given respectively by 


1 


S(x+) = i(1 — 2’)’ | oa) + 77 | P(x — t)f(t) dt 


— (| P(x — f(b) at) In (1 + »|, 


P ol 
S(x—) = —2i(1 — x)’ | ot) — m1 | P(x — bf(b dt 


— (| P(x — Of(b at) In (1 + |, 


where roots and logarithms of real positive quantities are to be interpreted in the con- 


ventional principal-value manner. It follows that 
al 
S(zt+) — S(a—) = (1 — 2’)' | oa) — (| P(x — df(t) at) In (1 + | 
\J -1 


and therefore, by a property of Cauchy integrals [5], 


S(z) = : (s — 2)"‘(1 — 8°)’ | 0) -- (| 
FW Ja} d 


P(s — Of( at) In (1 + »| ds + T(z), (6) 
1 
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where 7'(z) is entire. [7'\(z) could conceivably have isolated singularities at +1; however, 
the order of magnitude of the other terms in Eq. (6)—assuming that f(2) has no worse 
than a root singularity at these points— is such as to grow too slowly in the neighborhood 
of +1, so that no isolated singularity can exist.) Because the integral term in Eq. (6) 
tends to zero as z approaches infinity, 7'(z) must be the “entire part’’ of S(z)—..e., the 
positive power portion of the Laurent expansion of S(z) for the region outside the unit 
circle. Such a positive-power portion of a Laurent expansion can always be exhibited 
as a contour integral by means of Cauchy’s integral theorem, or alternatively as a residue 
at infinity; however, it is most convenient for our present purposes to use Eq. (4) directly. 

Consider next the quantity S(7+) + S(r—). Using Eq. (5) and again a property 
of Cauchy integrals [5] it is easily seen that 


1 


[ P(x — f(t) dt = h(2), (7) 


vz 


where 


h(x) = -Al-2x 2 | (s — x) ‘(1 — s°)' | a 
J -1 (8) 


al 


' -”= 
_ (| P(s — Of(b at) In (1 + »| ds + — rc}. 
= | Tr 

If now Eqs. (4), (7), and (8) are examined, it will be seen that although the unknown 
f(x) occurs on the right-hand side of Eq. (7), it does so only in the form of the (n + 2) 


constants 


al 


jo | f(t) dt, r=0,1,2,---,n4+1, 
« 1 

where n is the greater of the orders of the polynomials P and Q. Consequently, when 
Eq. (7) is solved for f(x) (see Sec. 3) it is only necessary to multiply the result by various 
powers of x and integrate between —1 and +1 in order to find linear equations relating 
these constants. Two such linear equations can always be obtained directly from Eq. 
(7) by multiplying both sides by (1 — «)' and setting z = +1 or —1. It may be verified 
that, with the exception of the integrals involving g(s), all integrals which occur in the 
process of setting up these various linear equations are easily expressed in terms of 
tabulated functions (in fact, nothing worse than Bessel functions). In this connection, 
the integrals tabulated in Sec. 5 are of value. 

In passing, an interesting alternative method of evaluating the c, may be mentioned. 
Suppose that the solution m,(x) were known for the equation 


| [Pit— a) Inj t—2|+Qt— DIA) dt = x. 
[f Eq. (1) is then multiplied by m,(x) and integrated with respect to x from —1 to +1, 


the result is 
1 al 


C= | Uf() dt = | m,(x)g(x) dx. 
J—1 -1 
3. Solution of Eq. (7). Taking Laplace transforms of Eq. (7) and inverting gives 


f(z) = a [ mae — an at, (9) 
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where M(x) is the inverse transform of [s°p,(s)]~', p,(s) being the transform of P(—.). 
In Eq. (8), A(x) has the apparent factor (1 — 2°) ! which would make for integration 
where h’ occurs. However, Eq. (7) shows that h(1) = Oand h(—1) 


difficulties in Eq. (9), 
must be 


finite (provided f(a 
cancelled by the remaining terms. This cancellation could be carried out explicitly if 
vr)? and setting wx +l and 


is integrable) and consequently the factor (1 


desired, by multiplying both sides of Eq. (7) by (1 
—1; if the first of these equations 1s multiplied by (4)(1 t+ x) (1 — 2°) ‘and the second 


and if the two resulting equations are added to Eq. (8), the 


by (3) (1 r) (1 
factor (1 r’)~? will be found to ecancel.* It will in fact be found that as x — 1, h(x) 
vanishes like (1 - this is easily seen to imply that f(x) behaves like (1 r) * as 


x — 1, and therefore similarly as « — —1 


2 Then a consideration of the 


4. First example. Set P? I, Q 0 [see Eq. (2)]. 
entire part of Eq. (4) shows that 7'(.x) is simply a constant: 


T(x) [ f(t) dt — | tf(t) dt. 


l 


These two constants, ¢ and Cc, , May be found by multiplying both sides of Kiq 7) by 


';) oe and setting +} and 1. The results are 


| [ q\s) ds 
) ( | 


: , ‘ I 
Tin 2 gs) 
| ‘ sq(s) ds 
C, | me 
©. (l — $s) 
Then using Iq. (9) gives 
d BP ae t Ss) “g(s) 
f(x) l ~ By | ds 
dr Ww J-1 8S—- x 
(10) 
af @ s°)'’" In (1 + 8s) ds l 
| (Cc - C,) 7 
W @— 4 T J 


From the fact that Eq 3) and (10) must be identical, the identity (11) can be deduced. 


5. Useful integrals. | 


the application of this method, a number of definite integrals 
ror convenience, the important ones are collected here. The first is 


continually occur; 


derived as just mentioned, and the others by straightfoward contour integration 

sf ts ; 

I n (1 d ea 
mt In 2 PG e ( — sin 11 

9 j 
| | | hn | ) ds 
: C(m + 1)r(r + 1 ; os F | 
2 ’ ; 3 fin2 + Yim + 1) — Wmt+re Ape 
C(m + r+ 2) 

*In carrving out the details of this discussion, it is advantageous to use an alternative form for 


T(z) found by use of a contour integral as previously mentioned, with the contour shrunk to envelop 


the strip 1, 1) 
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where because of the present rationality of m, r, the psi-functions are expressible in terms 
of elementary functions. for example, 
»l 7 1/2 
1+s 
(; In (1 + 8) ds = x(l — In 2), 
e _ == 


I. (; =) “In (1 +8) ds = —x(1 + In 2), 


1 
. 


*! sF(s) ds ls) ds 
| Me) a = | F(s)ds + x Fs) ds ; 
1 §$—- 72@ J-1 4 3 | 
| ( = 3° 1/2 Is vl 1 iad 3°) 1/2 Is 
l on iti | ( s) “ds - 
1 g=— fs 1 a= 2 
on (2 oe ee s)'’* ds pt ofl — g) *”* de er (1 +s '7(1 — 8)'" ds 
| paoe., _pwsermeae 
P , Ss MH 1 o=—= 2 e 1 ee | 
d [ F(s) ds 7 Ful) F(—1) 4 [ F'(s) ds 
dz J-; 1 — 2 l—@z l+e2 Jan @—-2z’ 


I'(n + 1/2)1,(—a) 


e“(1 — #)""" dt a 
I x '’*(—1/2a) 


(J, = modified Bessel function), 


| i sin : idt 3 [cosh es I,(a)], 
F a 


1 " 1 
. ; : : - . ta : 
| t"e* F(b dt - nz | e*' F(t) dt. 
J-1 da . 1 
6. Second example. Consider 


1 
. 


| (+ ke - Od] |x — ¢| fade =1, (12) 


where k is a parameter. It is of interest to compare the solution of Eq. (12) with that 
obtained by setting k = 0 and using Eq. (3); in particular, one may ask how large k 
must be before its effect becomes important. The function 7'(x) is found from Eq. (4) as 


T(x) = eo(—1 + 1/2k) — kx] + ¢,[(—1 + &) — kr] + 1/2ko, 
and two relations between these ec, may be found by the method of See. 4: 
co(In 2) + ke,(1 — m2) = —-1, 
keo(1/2 — In 2) — 2c, + ke, = 0. 


These two equations may be used to simplify the expression for h(x), which becomes 
(had ¢ sew I : 
h(x) [eo(—k In 2)] + (5 — sin ' rfc + kx) — ke,]}. 
7 T \“ 
linally, Eq. (9) gives 
Boa d f' a 
f(z) = = | exp [A(t — x)|h'(0) de 
_ (13) 


1 
. 


—h'(x) + kh(x) + | exp [a(t — x)]h(0) dt. 


vr 
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Integrating Eq. (13) between —1 and +1 gives another relation between the c, 
¢o(1 — In 2)/,(k) — e,1.(k) = 0 


so that, solving for c, and ¢, , 


[,(k) (1 — In 2)/,(k) 
oO =- ———— ‘ ¢°¢=->z ree 
(In 2)Io(k) + kI,(k) (In 2)I,(k) + kI,(k) 
and the final solution may now be written down by substituting for h(x) and the ¢; in 
Eq. (13). The special case of small / is of interest; to the first order, 
; ] ( l ) £ 
f(x) = — - “773 kl — —— “ot ie 
: wr In 2(1 — 2°)'” * In 2/ (1 — a*)'” 
where the first term is the well-known solution of Eq. (12) for k = 0. The correction term 


bears a ratio to the first term of approximately (3kx), which is of the expected magnitude. 

To conclude, it may be remarked that to write out a formal solution is estimated to 
require about 2n*’* hours of labor for that case in which P and Q are nth order poly- 
nomials. Of course, once such a solution has been tabulated, it may be used for a number 
of kernels or parameters. Note also, as previously mentioned, the possibility of curve- 
fitting. In particular, if a large number of terms are to be used, it may be worthwhile 
to set P = 1 and represent the kernel behavior by [In | x — ¢! +Q] alone; the Q terms 


could then of course be absorbed directly into g(x) and Eq. (3) used. 
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PRIORITY ASSIGNMENT ON A WAITING LINE* 


By 8. A. DRESSIN! anp E. REICH? (Santa Monica, California) 


Interest in queueing statistics when a priority system is involved is now widespread, 
especially in connection with communications applications. In this paper we consider 
a model described in the usual customer-counter terminology as follows. 

Input. There are r classes of customers, a particular class, labelled by the integer p, 
(p = 1,2, ... ,7), being made up of customers of ‘‘priority” p. (I’ollowing an unfortunate, 
but standard, practice, priority p, will be said to be higher than priority p, if and only 
if p, < p2 .) Customers of priority p arrive as a Poisson process with interarrival-period 
distribution \,e ‘dx , x > 0. Customers of different priorities arrive independently. 

Queue-discipline. Once a customer’s service has begun, it is permitted to proceed to 
completion. If the counter becomes empty, and there is at least one customer waiting, 
*Received May 7, 1956. 

Major, U. 8. M. C., U. 8. Naval Postgraduate School, Monterey, Calif. 
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then a customer of the highest priority present on queue is admitted to the counter, and 
his service starts. Customers of equal priority are served first come, first served. 

Service mechanism. The service period has the distribution we “dz, x > 0. Service 
periods are mutually independent, and independent of input history and priority. 

The above model will be recognized to be a specialization of one considered by 
Cobham [1]. It is reasonable in some applications, and has the mathematical virtue of 
yielding simpler results than more general queue disciplines. This paper extends previous 
work in that we compute conditional and unconditional distributions, rather than only 
means. 

In what follows, we put 


> 
o, = a & o,=40. 


t=1 


‘The system will be said to be in state n if the total number of customers at the counter 


and on queue is n(n = 0, 1, 2, ---). At the same time, the system will be said to be in 
state (n; p), n = 1, 2, --- , if there is a customer of arbitrary priority at the counter, 
and there are exactly n — 1 customers, of priority p or higher, on queue. State (0; p) 


will be identified with state 0. 

System behavior with time. It is clear that the value of n at time ¢ is the same as 
if we were dealing with the classical case r = 1. The reason for this is that the particular 
queue discipline under consideration can be transformed to first come, first served, by 
merely reordering those customers on queué, in a manner not depending on their service 
times. If we continue this reasoning, applying it to the class of cutomers of priority 
p or higher, we see that the following is true: For every fixed p, the vector A = [k, (m; p)] 
is a Markoffin function of time. 

Let P,(t) = probability that the system is in state n at time ¢, P,,.,,)() = probability 
that the system is in state (n; p) at time ¢. In view of what was said above, the differential 
equations for P,,(2) are the well-known ones [4], 


0.P,(t) = —oP(t) + uPi(d, (1) 
0.P,(t) = oP,-(0 — (o + w)P,(0) + uP, (0, n=i1,3,--++. (2) 

Trivially, 
P.(f = P£O. (3) 


Let h > 0 be a small quantity. Neglecting probabilities of order h*, we can say that 
the system is in state (1; p) at time ¢ = ¢, if and only if at least one of the following 
mutually exclusive events has occurred. 


(i) State at / = ¢, — h was (1; p), no service ended since 4, — h, no customers of 
priority p or higher arrived since 4, — A, or 

(ii) the value of A at time f, — h was [k, (1; p)], & > 1, a service period ended since 
f, — h, no customers of priority p or higher arrived since f, — h, or 

(iii) state at 4, — h was (0; p), a customer of arbitrary priority arrived since f, — h, or 

(iv) state at 4, — h was (2; p), a service period ended since ¢, — h, no customers of 


priority p or higher arrived since f, — h. Now, within terms of order h’, 
Prob (i) = (1 — ph — o AP a.) (to — A), 
Prob (ii) = wh[P,,.,(t — h) — P(t — Ad], 
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Prob (iii) ohP(t - h), 


Prob (iv) BAP (2: (t — h). 
Hence 
6.F (Z) oP, (t) — uP) —- Gok cast { uP o.»)(b). (4) 
Similarly, 
OF tavnsit GF, t) — (o, + wPiip() + wPinsi:y) (0; n= 2,3, °°: (5) 


Note that Eqs. (3)-(5) involve P,(4), P,(4), as determined by Eqs. (1), (2). Ifo < yz 


then the “equilibrium” solutions P, , P , of (1)-(5), for which 


Pe. win RK. 
are 


P (# *\(2) , k=0,1,2,°°:, (6) 


a mu 
and, therefore 

Li o L 0,)o ‘og j ny 
P P 4 ele (22) , k=1,2,---. (7) 

M MO, bu 

Ifo > uw > o,, the appropriate values are 
P; 0, (k OG, Bs me 
P a) Po. (# 0 )( 2 (} i ee © 


Conditional queueing times. We now calculate the probability f;,.,)(7)d7T that a 
customer C of priority p will have to wait on queue for a duration of between 7’ and 
T + dT, conditional to the assumption that the system is in state (n; p) the instant 
before C’s arrival. It is clear that the probability in question would not be altered if it 
were further conditioned relative to some aspect of the queue’s history prior to C’s 
arrival. 

Of course 

T’) (7 0), p Pe are (8) 


Also, 


T n fold convolution of we *’ ane. Sm 1)!, (9) 


By employing a simple device we can immediately write down f(;;,)(7'), p > 1. 
If a customer C of priority p finds, on arrival, that the only customer delaying him is a 
customer at the counter, then C’s queueing time 7’ depends on the additional service 
time of the customer at the counter, plus the service times of the customers of priority 
higher than p who arrive while the customer at the counter is served, etc. 7 evidently 
has the same distribution as the length of the busy period of the classical single-priority 


queue with input and service parameters ¢,_, and y, respectively, that is, [3, 4], if 
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> 


 < H Saw (T) = @/e, as oll exp [—(¢,-. + wT VI, (20)) p' *T'), 
T>0, p=2,3,-:-, 

where /,, is the modified Bessel function of order n. If the system is in state (n; p) n > 1, 
p > 1, the instant before C’s arrival, say at ¢ = 0, then 7’ is unaffected if we reduce the 
priorities of those customers of priority higher than p, who are on queue at / = 0, to p. 
Assume this reduction is carried out, Let C* denote that customer of priority p who is 
directly in front of C, and Jet 7* denote the remaining queueing time of C*, from ¢ = 0 
on. 7* has the distribution f,,,-,,,)(a)dz. At t = 7*, C*’s service starts, and C is first 
on line in front of the counter. Hence the variable 7’ — 7* has distribution f,,.,)(x)dz, 
and is independent of 7*. Consequently, f,,.,)(«) is the convolution of f,,-,.,)(@) and 
fa:p) (x); that is, 
f (7) n fold convolution of f(,.,)(7) 

n(u/o,-1)""T * exp [—(e,-1 + wT IL, (20,707), (10) 
3 


T > 0, n ee p = 2,3, °°: 


Unconditional queueing times. Suppose test-customer of priority p enters the 
system at time /. The probability that he will have to wait on queue for a period between 
T and T + dT is 


FAT; 0 AT = YS Pay (Of a» (T) aT. 
k=-0 
In “equilibirum’’, we have, by (7)-(10), if o < u, 
fe") = 1 — (6/p) + 2(o/n)(u — o,)Q™', (11) 


W here 


1/2 1/2\1/2 ‘ 2 ? ) 
8+ 0,-; tut 2o, mu) (8 + op-1 + — 2o,-e”) 


i) S+o Lu — 2a, + 


« 


a = 0. 


\n analogous expression is obtained when ¢ > yu > "a 
Note that Eq. (11) gives, for p = 1, 


m u-d . 
Prob {T < 7} l — (o/pe ae zg > 0. 
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THE RELATION BETWEEN THE FLOW OF NON-NEWTONIAN FLUIDS AND 
TURBULENT NEWTONIAN FLUIDS* 


By R. S. RIVLIN (Brown University) 


It is well known [1]** that if a Newtonian fluid such as water flows through a tube 
of circular cross-section under conditions for which the critical Reynolds number is 
exceeded, the relation between pressure gradient and volume rate of flow is non-linear. 
Further [2], the radial distribution of velocity (averaged over regions large compared 
with the eddy dimensions) departs markedly from the parabolic distribution which 
would be expected for a Newtonian fluid. These results suggest that some analogy may 
exist between the turbulent Newtonian fluid and the laminar flow of a non-Newtonian 
fluid. It is the object of this note to pursue this line of thought and to indicate some of 
the phenomena which might, in accordance with it, be observed in turbulent. fluids. 

If we assume that in a fluid, which is isotropic in its state of rest, the stress com- 
ponents in the rectangular Cartesian coordinate system x; may be expressed as poly- 


nomials in the gradients of velocity, acceleration, second acceleration, --- , (n — 1)th 
acceleration, then the stress matrix T (= || ¢;; ||) may be expressed as a matrix poly- 
nomial in n kinematic matrices A, , A, , --- , A, , defined by [3, 4] 
Ov; Ov, 
A, = || Ai; od | e-em sar on 
OX; OX; 
and 
- r 1 ? 1 
OA; aA ' a on 
A, = || A fi 9, + AS? = a. = 
ot Ox Ox Ox, (1) 
(y = 2.3, , n) 


where v; is the velocity of the fluid. The coefficients in this matrix polynomial are poly- 
nomial invariants, under orthogonal transformations, of A, , A, , --: , A,. If A, = 0 
for r > 2, then T may be expressed as a matrix polynomial in the two matrices A, and 
A, . It has been shown [5] that it may then be expressed in the form 

T = —plI + 0,A, + aA, + a;A; + a,A; + a;(A,A, + A,A,) 

+ a,(AiA, + A.A}) + a;(A,A; + A°A,) + a.(ATA; + A;A)), (2) 
where p is arbitrary if the velocity field is specified, I is the unit matrix and the a’s are 
polynomials in ¢rA{ , trA? , trA{ , trA} , trA,A, , trA{A, , trA,A? , trA{A? . In deriving this 
result it is assumed that the fluid is incompressible. 

It can be shov n [6] that the assumption that the stress components at a point of the 
fluid can be expressed as polynomials in the gradients of the velocity, acceleration, etc. 
at that point is, under rather general conditions, equivalent to an assumption that the 
stress in an element of the fluid at any instant of time depends on the velocity gradients 


in the element at that instant of time and at previous instants. 





*Received June 25, 1956. The results presented in this note were obtained in the course of research 
sponsored by the Office of Ordnance Research, U. 8S. Army, under Contract No. DA-19-020-3487. 
**Numbers in square brackets refer to the Bibliography at the end of the paper. 
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We shall now consider the results obtained on the basis of this theory for the flow 
of a fluid along a straight pipe of non-circular cross-section. If we assume that each 
particle of the fluid moves in a rectilinear path parallel to the length of the pipe (in the 
x3-direction say), so that 


i =e, = 0 and vv = 0,(2;, 2), (3) 


then it is easily seen, from (1), that A, = 0 (r > 2), so that Eq. (2) is applicable. We 
now substitute the expressions (3) for the velocity components in (2) and then substitute 
the resulting expressions for the stress components in the equations of motion with 
zero body forces. The three equations so obtained may be regarded as differential equa- 
tions for the determination of p and v, , subject to the boundary condition that v, = 0 
on the wall of the pipe. It has been shown by Ericksen [7] that these equations are 
generally incompatible unless the pipe is circular. (Exceptions occur for particular 
forms of the expression (2) for the stress; for example, the form T = —pI + 7A, with 
n constant, which represents the Newtonian fluid, provides such a case.) This incom- 
patibility implies that for a non-Newtonian fluid to which Eq. (2) applies, no state of 
rectilinear flow along a non-circular pipe is in general possible, unless body forces are 
applied to the fluid. Green and Rivlin [8] have considered the particular case of a nearly 
Newtonian fluid for which the stress is given by 


T = —pl + nA, + ec + trAA; , (4) 


where 7 and ¢ are constants and ¢ is small enough so that any departure from the flow 
field predicted for the Newtonian fluid, obtained by taking « = 0 in (4), may be treated 
as a first order perturbation. They have found, for the case of a pipe with elliptical 
cross-section, a flow field in which secondary flows in the cross-sectional planes are 
superimposed on the rectilinear flow. Further calculations with somewhat similar results 
have been carried out by Langlois and Rivlin [9] using a more general form than (4) 
for the constitutive equation of the nearly Newtonian fluid. 

It has been reported [10] that when a Newtonian fluid flows down a straight pipe of 
non-circular cross-section, under conditions for which the fluid has become fully turbu- 
lent, the flow is no longer rectilinear, but a secondary flow exists in the cross-sectional 
planes of a type similar to that calculated by Green and Rivlin and by Langlois and 
Rivlin. This fact further suggests that the turbulent Newtonian liquid may, for certain 
purposes, be regarded as a non-Newtonian fluid. 

On the basis of the theory for the flow of non-Newtonian fluids described above, 
we can derive the existence of other phenomena. If the fluid is in a state of flow with 
uniform velocity gradient K, described by 


o, = Kz, , vy, =v, = 0, (5) 


then, from (5), (1) and (2), A, = 0 (r > 2) and the stress components /,;; at each point 
of the fluid are given by [11] 


ti, = —p+a,K’, 

tee = —p + K*[(2a, + as) + 4(ay + a%)K* + 8asK"], (6) 
tio = Kla, + 2a;K? + 4a;K"‘], 

ts, = —p, beg = ty, = 0. 
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Unlike the situation that exists in a Newtonian fluid, the normal components of stress 
are not equivalent to a hydrostatic pressure. This results [11] in a number of phenomena 
which have been observed in non-Newtonian fluids and may be called normal stress 
effects. 

For example, it was remarked by Garner and Nissan [12] that, for certain fluids, if 
a cylindrical rod is rotated in the fluid, the fluid rises up the rod. For a Newtonian 
fluid, the fluid surface would remain horizontal under conditions for which the centrifugal 
forces are negligible and be depressed slightly at the stirrer if centrifugal forces were not 
negligible. The rise at the stirrer which occurs with the non-Newtonian fluids implies 
that a radial distribution of normal surface thrusts would have to be exerted on the 
free surface of the non-Newtonian fluid in order to maintain in it a velocity field in 
which each particle moves with uniform velocity in a horizontal circular path about 
the axis of rotation. 

Again, Weissenberg [13] contained certain non-Newtonian fluids between two coaxial 
the clearance between the bases of which was small compared 


] 


¢evlindrical containers, 
with their diameters. The outer cylinder was rotated at a constant angular velocity, 
while the inner one, to which were attached a number of vertical tubes communicating 
with the fluid through small holes, was held stationary. He found that the fluid rose in 
the tubes to an extent depending on their distance from the axis of rotation, the rise 
being greatest at the axis and decreasing with increasing radial distance. In the region 
between the dises, each particle of the fluid moves substantially in a horizontal circular 
path about the axis of rotation with a uniform velocity depending on the radial and 
vertical positions of the particle. The rise of fluid in the tubes implies that, in order to 
maintain this state of flow in the non-Newtonian fluids, surface tractions must be exerted 
on the plane surfaces of the fluid mass directed normally into the fluid. For a Newtonian 
fluid, of course, only shearing tractions are required in order to maintain a similar flow*. 

It can be shown that the terms in the expressions for the stress components, which 
give rise to the secondary flow in the problem of the flow of a non-Newtonian fluid 
down a non-circular pipe mentioned above, also give rise to normal stress effects of the 
types observed by Garner and Nissan and by Weissenberg and this raises the possibility 
that similar normal stress effects may also be observed in turbulent Newtonian fluids 
in analogous situations 

It has been noted that the fluids in which normal stress effects have been observed 
are visco-elastic and this suggests that turbulent Newtonian fluids may similarly show 
some elastic effects under suitable experimental conditions. It may be remarked that in 
solutions, which exhibit both normal stress and visco-elastic effects, the 


high-polym 
origin of the phenomena probably lies in the preferential orientation and extension ot 
the dissolved high-polymer molecules when the liquid is sheared. The eddies in a turbulent 
Newtonian fluid will presumably undergo preferential orientation when the turbulent 
fluid is sheared providing a possible mechanism for the effects in the turbulent fluid 


Also, the eddies in the sheared volume elements of the turbulent fluid would presumably 


*This is not strictly correct if the effect of centrifugal forces is taken into account, but centrifugal 
forces would cause the fluid to rise in the vertical tubes increasingly with increasing radial distance 
from the axis of rotation. Furthermore, the effect observed with non-Newtonian fluids may be quite 


pronounced under conditions of negligible centrifugal force. 
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be distributed in an anisotropic manner and this would also provide a possible mechanism 
for the effects discussed. 

It will, of course, be appreciated that the above remarks are of a qualitative nature 
and are presented as tentative suggestions rather than as positive predictions. Certainly 
it appears quite likely that a phenomenological theory of the type considered, in which 
the stress in an element of the turbulent fluid (large compared with the eddy dimensions) 
is supposed to depend only on the kinematic variables in that element or on the velocity- 
gradient history of that element, will be entirely adequate as a complete description of 
the flow properties of the turbulent fluid, since eddies can diffuse from one point of the 
fluid to another. 
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Infinite sequences and series. By Konrad Knopp. Translated by Frederick Bagemihl. 
Dover Publications, Inc., New York, 1956. 186 pp. $3.50. 


In this little volume, about one fourth as long as the author’s ‘“Theory and Application of Infinite 
Series,’’ Knopp presents an introduction to the theory of infinite series. With consummate skill he leads 
the reader through the elementary notions and then well into the more subtle developments. The chapter 
headings are as follows: Introduction and prerequisites, Sequences and series, The main tests for infinite 
series-operating with convergent series, Power series, Development of the theory of convergence, Ex- 
pansion of the elementary functions, Numerical and closed evaluation of series. With a few exceptions, 
the results presented are included in the treatise mentioned above. In particular, the following are new: 
Jehle’s extension (1950) of the Weierstrass convergence test (concerning the asymptotic form of an4:/a, 


for a series > a,); Landau’s proof (1920) of Hardy’s theorem on multiplication of series. 
WILFRED Kaplan 


Selecta Hermann Weyl. Published on the occasion of his seventieth birthday by the 
Federal Institute of Technology in Ziirich and the Institute for Advanced Study in 
Princeton. Birkhauser Verlag, Basel and Stuttgart, 1956. 592 pp. $11.45. 

To celebrate the seventieth birthday of Hermann Wey], on the ninth of November 1955, the Federal 

Institute of Technology in Ziirich and the Institute for Advanced Study in Princeton together published 

this commemorative volume. These two institutions, to whom Hermann Wey! dedicated the greatest 


part of his scientific activity, wished to express their appreciation of his work in research and teaching 


and their admiration of his mathematical genius. 

This collection of articles, written during the period from 1910 to 1952, is a cross-section of Weyl’s 
mathematical work, and reflects the advances in almost all fields of mathematics during this interval 
of time. The text of the papers included in this collection has been essentially unaltered, save for minor 
changes in the style of the original bibliographical references and the correction of obvious misprints. 
Many of the papers reappear here with valuable annotatians and cross references (dated 1955 )by Hermann 
Weyl himself. The first paper is the well known 1910 Mathematische Annalen article dealing with 
“singular” ordinary differential equations d/ds (p(s)du/ds) — q(s)u(s) = 0, in which Weyl introduced 
his “limit point’’ and “‘limit circle” criteria. The second paper is the 1915 Rendiconti di Palermo article 
on the distribution of the eigenvalues of an arbitrary elastic body. This paper contains the author’s 
only major alteration of the original text, in connection with the proof of the existence of a solution of a 
certain non-homogeneous integral equation whose corresponding homogeneous integral equation may 
possess non-trivial solutions. The last paper in the book, Mathematische Zeitschrift 1952, deals with 
boundary value problems of the same nature (but in a more general setting) as those considered in the 
second paper just mentioned. There are nineteen papers in all, arranged chronologically. Special mention 
will be made here only of the 1916 Mathematische Annalen paper containing Weyl’s generalization of 
Kronecker’s approximation theorem; the 1916 Vierteljahrschrift der Naturforschenden Gesellschaft 
in Ziirich paper on the determination of a closed convex surface by its line element; the 1919 Annalen 
der Physik paper on the propagation of electromagnetic waves; the 1927 Mathematische Annalen paper 
on integral equations and almost periodic functions; the 1935 American Journal of Mathematics paper 
(with R. Brauer) on spinors in n dimensions; the 1940 Duke Mathematical Journal paper on the method 
of orthogonal projection in potential theory; and the 1942 Annals of Mathematics paper on the differ- 
ential equations of the simplest boundary layer problems. This incomplete list of the contents is meant 
to convey an idea of the rich source of inspiration embodied in this single volume. The book concludes 
with a complete list of the publications of Hermann Weyl. 

The appearance of a volume of this kind is always a welcome event in the scientific world. It is to 
be hoped that this beautifully appointed Birkhauser Verlag volume is but a harbinger of the publication 


of the collected works of Hermann Weyl. 
J. B. Diaz 




















